如何解决在Python中使用Runge-Kutta求解耦合微分方程组的系统
import numpy as np
import matplotlib.pyplot as plt
import numba
import time
start_time = time.clock()
@numba.jit()
# A sample differential equation "dy / dx = (x - y**2)/2"
def dydx(x,y):
return ((x - y**2)/2)
# Finds value of y for a given x using step size h
# and initial value y0 at x0.
def rungeKutta(x0,y0,x,h):
# Count number of iterations using step size or
# step height h
n = (int)((x - x0)/h)
# Iterate for number of iterations
y = y0
for i in range(1,n + 1):
"Apply Runge Kutta Formulas to find next value of y"
k1 = h * dydx(x0,y)
k2 = h * dydx(x0 + 0.5 * h,y + 0.5 * k1)
k3 = h * dydx(x0 + 0.5 * h,y + 0.5 * k2)
k4 = h * dydx(x0 + h,y + k3)
# Update next value of y
y = y + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4)
# Update next value of x
x0 = x0 + h
return y
def dplot(start,end,steps):
Y=list()
for x in np.linspace(start,steps):
Y.append(rungeKutta(x0,y,h))
plt.plot(np.linspace(start,steps),Y)
print("Execution time:",time.clock() - start_time,"seconds")
plt.show()
start,end = 0,10
steps = end* 100
x0 = 0
y = 1
h = 0.002
dplot(start,steps)
dydx= (x - y**2)/2
dydt= (x - y**2)/2
dxdt= x*3 + 3y
如何在上述代码中将这两个实现为耦合微分方程组? 具有n个耦合微分方程组的系统是否还有更通用的方法?
解决方法
在其他人的帮助下,我做到了:
import numpy as np
from math import sqrt
import matplotlib.pyplot as plt
import numba
import time
start_time = time.clock()
a=1
b=1
c=1
d=1
# Equations:
@numba.jit()
#du/dt=V(u,t)
def V(u,t):
x,y,vx,vy = u
return np.array([vy,a*x+b*y,c*x+d*y])
def rk4(f,u0,t0,tf,n):
t = np.linspace(t0,n+1)
u = np.array((n+1)*[u0])
h = t[1]-t[0]
for i in range(n):
k1 = h * f(u[i],t[i])
k2 = h * f(u[i] + 0.5 * k1,t[i] + 0.5*h)
k3 = h * f(u[i] + 0.5 * k2,t[i] + 0.5*h)
k4 = h * f(u[i] + k3,t[i] + h)
u[i+1] = u[i] + (k1 + 2*(k2 + k3 ) + k4) / 6
return u,t
u,t = rk4(V,np.array([1.,0.,1.,0.]),10.,100000)
x,vy = u.T
# plt.plot(t,x,t,y)
plt.semilogy(t,y)
plt.grid('on')
print("Execution time:",time.clock() - start_time,"seconds")
plt.show()
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。