Я написал класс с целью решения системы дифференциальных уравнений (данной в форме numpy.array), для решения нелинейной системы я использую scipy.optimize.fsolve, используя пример, найденный здесь в одном посте, метод отлично работает с одним уравнением, но потерпит неудачу, если я попытаюсь использовать для системы дифференциальных уравнений! Я написал минимальный, полный и проверяемый пример, чтобы вы могли проверить и глубоко понять, как работает класс!
import numpy as np
from scipy.optimize import fsolve , newton_krylov
import matplotlib.pyplot as plt
class ImpRK4 :
def __init__(self, fun , t0, tf, dt , y0):
self.func = fun
self.t0=t0
self.tf=tf
self.dt=dt
self.u0=y0
self.n = round((tf-t0)/dt)
self.time = np.linspace(self.t0, self.tf, self.n+1 )
self.u = np.array([self.u0 for i in range(self.n+1) ])
def f(self,ti,ui):
return np.array([functions(ti,ui) for functions in self.func])
def solve(self):
for i in range(len(self.time)-1):
def equations(variable):
k1,k2 = variable
f1 = -k1 + self.f(self.time[i]+ (0.5+np.sqrt(3)/6)* self.dt , self.u[i]+0.25*self.dt* k1+ (0.25+ np.sqrt(3)/6)*self.dt*k2)
f2 = -k2 + self.f(self.time[i]+ (0.5-np.sqrt(3)/6)* self.dt , self.u[i]+(0.25-np.sqrt(3)/6)*self.dt *k1 + 0.25*self.dt* k2)
return np.array([f1,f2]).ravel() #.reshape(2,)
k1 , k2 = fsolve(equations,(2,2)) #(self.u[i],self.u[i]))
self.u[i+1] = self.u[i] + self.dt/2* (k1 + k2)
plt.plot(self.time,self.u)
plt.show()
def main():
func00 = lambda t,u : -10*(t-1)*u[0]
func01 = lambda t,u : u[1]
func02 = lambda t,u : (1-u[0]**2)*u[1] - u[0]
func0x = np.array([func00])
func0 = np.array([func01,func02])
t0 = 0.
tf = 2.
u0 = y01
dt = 0.008
y01 = np.array([1.,1.])
diffeq = ImpRK4(func0,t0,tf,dt,y01)
#y0 = np.array([np.exp(-5)])
#diffeq.solve()
#diffeq = ImpRK4(func0x,t0,tf,dt,y0) ## with single equations works
diffeq.solve()
if __name__ == '__main__':
main()
EDIT Нет, извините, но это не то, что я искал... в основном, когда у меня есть система уравнений, я должен получить K1 и K2 той же размерности self.u[i]
