Как изогнуть три кривые экспоненциального затухания с разными параметрами? - PullRequest
Новая
       18

Как изогнуть три кривые экспоненциального затухания с разными параметрами?

0 голосов
/ 09 июля 2019

Я не слишком знаком с Python, но пытался изогнуть три разных энергетических кривых классического гармонического осциллятора. Я нашел один метод, который не самый эффективный, если я хочу подгонять кривую больше. 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy import optimize
import sympy as sy

dU=np.zeros(2)                
def oscil(U,t):            
    dU[0]=U[1]
    dU[1]=-U[0]-b*U[1]
    return dU

ic = [0,1]                  #initial condition
t = np.linspace(0,100,1000) #time

#4 different graphs with different beta, B, values; B=[0,0.1,0.5,1]
b=0
soln1 = odeint(oscil, ic, t)

b=0.1
soln2 = odeint(oscil, ic, t)

b=0.5
soln3 = odeint(oscil, ic, t)

b=1
soln4 = odeint(oscil, ic, t)

#Create 4 subplbots
fig, ( (plt1,plt2), (plt3,plt4) ) = plt.subplots(2,2)
fig.subplots_adjust(wspace=0.4, hspace=0.8)
plt1.plot(t, soln1[:,0], label='position')
plt1.plot(t, soln1[:,1], label='velocity')
plt1.set_title('Undamped Harmonic Oscillator')
plt1.set_xlim(-1,4*np.pi)
plt1.set_xlabel('Time')
plt2.plot(t, soln2[:,0], label='position')
plt2.plot(t, soln2[:,1], label='velocity')
plt2.set_title('Damped HO for B=0.1')
plt2.set_xlabel('Time')
plt3.plot(t, soln3[:,0], label='position')
plt3.plot(t, soln3[:,1], label='velocity')
plt3.set_title('Damped HO for B=0.5')
plt3.set_xlim(-1,50)
plt3.set_xlabel('Time')
plt4.plot(t, soln4[:,0], label='position')
plt4.plot(t, soln4[:,1], label='velocity')
plt4.set_title('Damped HO for B=1')
plt4.set_xlim(-1,20)
plt4.set_xlabel('Time')
plt.title('HMO with Driving Force')
plt.legend(loc='best')
plt.show()

# Energy = x^2 + dxdt^2
Energy1 = soln1[:,0]**2+soln1[:,1]**2
Energy2 = soln2[:,0]**2+soln2[:,1]**2
Energy3 = soln3[:,0]**2+soln3[:,1]**2
Energy4 = soln4[:,0]**2+soln4[:,1]**2

# curve fit
x=t 
y2=Energy2 
y3=Energy3 
y4=Energy4 

def mod2(x, xshift, steepness, yshift):  ###estimated function
    return xshift*np.exp(-steepness*x) + yshift

init_vals2 = [1.11913923, 0.11376818, 0.0826447 ]  
init_vals3 = [1.21244276, 0.57442609, 0.07596321 ]  
init_vals4 = [1.04292458, 0.90770071, 0.06833576]       
                
fitParams2, fitCovariances = optimize.curve_fit(mod2, x, y2, p0=init_vals2)
fitParams3, fitCovariances = optimize.curve_fit(mod2, x, y3, p0=init_vals3)
fitParams4, fitCovariances = optimize.curve_fit(mod2, x, y4, p0=init_vals4)
plt.plot(t,y2, t,y3, t,y4)
plt.plot(t,mod2(x,*fitParams2), 'k' )
plt.plot(t,mod2(x,*fitParams3), 'k' )
plt.plot(t,mod2(x,*fitParams4), 'k' )
plt.title('Energy')
plt.xlabel('t')
plt.show()

Я бы хотел что-то вроде: 

y = [Energy2, Energy3, Energy4]
mod=[mod2, mod2, mod2]
p0= [init_vals2,init_vals3,init_vals4]


for i in range(0,2):
    fitParams, fitCovariances= optimize.curve_fit(mod[i], x, y[i], p0[i])

for i in range (0,3):
    plt.plot(t,y[i]) 

for i in range (0,2):
    plt.plot(t, mod[i+1](x, *fitParams), 'k')
plt.show()

Однако, он может соответствовать только одной кривой. Я не могу найти какие-либо решения в другом месте. Спасибо!

...