Как интегрировать связанные дифференциальные уравнения?

Question

У меня есть система уравнений, которую я пытался заставить Python решить и построить, но сюжет получился неправильным.Это мой код:

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


#function that returns dx/dt and dy/dt
def func(z,t):
    for r in range(-10,10):
                beta=2
                gamma=0.8

                c = z[0]
                tau = z[1]
                dcdt = r*c+c**2-c**3-beta*c*tau**2
                dtaudt = -gamma*tau+0.5*beta*c*tau
    return [dcdt,dtaudt]

#inital conditions
z0 = [2,0]

#time points
t = np.linspace(0,24,100)

#solve ODE
z = odeint(func,z0,t)

#seperating answers out

c = z[:,0]
tau = z[:,1]

print(z)


#plot results
plt.plot(t,c,'r-')
plt.plot(t,tau,'b--')
plt.legend(['c(t)','tau(t)'])
plt.show()

Позвольте мне объяснить.Я изучаю вдвойне диффузную конвекцию.Я не хочу делать какие-либо предположения относительно значения r, но бета и гамма положительны.Поэтому я хотел бы присвоить им значения, но не r.Это сюжет, который я получаю, и из понимания проблемы, что график не верен.График тау не должен зависеть от 0, а график c должен делать больше.Я плохо знаком с Python и прохожу курсы, но очень хочу понять, что я сделал неправильно, поэтому помощь на простом языке была бы признательна.

Answer 1

Я вижу 2 проблемы в вашей функции, которые вы должны проверить.

    for r in range(-10,10):

Здесь вы выполняете цикл for, просто переоценивая dcdt и dtaudt.В результате выходное значение совпадает с простым вычислением r = 9 (последнее значение в цикле)

    dtaudt = -gamma*tau+0.5*beta*c*tau

Здесь у вас есть dtaudt = tau*(beta*c/2. -gamma).Ваш выбор tau[0]=0 подразумевает, что тау останется 0.

Попробуйте это:

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

r = 1
beta=2
gamma=0.8

#function that returns dx/dt and dy/dt
def func(z,t):
    c = z[0]
    tau = z[1]
    dcdt = r*c+c**2-c**3-beta*c*tau**2
    dtaudt = -gamma*tau+0.5*beta*c*tau
    print(dtaudt)
    return [dcdt,dtaudt]

#inital conditions
z0 = [2,0.2] #tau[0] =!0.0

#time points
t = np.linspace(0,24,100)

#solve ODE
z = odeint(func,z0,t)

#seperating answers out

c = z[:,0]
tau = z[:,1]

#plot results
plt.plot(t,c,'r-')
plt.plot(t,tau,'b--')
plt.legend(['c(t)','tau(t)'])
plt.show()