[已解决]微分方程的龙格-库塔法Python编程请教

kevinheros

有个微分方程要用龙格-库塔法求解，但我编程求解后的图和实际的图不一样，是否有了解的鱼油看看问题所在？多谢！！！
1. 原微分方程：

2.边界条件：（其中，h=8mm，理解为当x=0时，y=8）

3.别人用Matlab编程得到的曲线图：

4.我用Python编写的程序及对应的曲线图：

1. import numpy as np
   
2. import matplotlib.pyplot as plt
   
3. 
   
4. # 龙格-库塔法函数定义
   
5. def runge_kutta(y, x, dx, f):
   
6.     """ y is the initial value for y
   
7.         x is the initial value for x
   
8.         dx is the time step in x
   
9.         f is derivative of function y(t)
   
10.     """
    
11.     k1 = dx * f(y, t)
    
12.     k2 = dx * f(y + 0.5 * k1, x + 0.5 * dx)
    
13.     k3 = dx * f(y + 0.5 * k2, x + 0.5 * dx)
    
14.     k4 = dx * f(y + k3, x + dx)
    
15.     return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6.
    
16. 
    
17. if __name__=='__main__':
    
18.     t = 0.0
    
19.     y = 8.0 # 初始值为(0,8)
    
20.     dt = 0.000001 #步长
    
21.     ys, ts = [], []
    
22. 
    
23.     # 微分方程函数定义
    
24.     def func(y, t):
    
25.         v1 = np.sqrt(t**2 + y**2)
    
26.         v2 = 1000 * t / v1
    
27.         v3 = np.sqrt((v2 - t)**2 + (1000 - y)**2)
    
28.         return (v1*(v2-t)-1.49*t*v3)/(1.49*y*v3-(1000-y)*v1)
    
29.    
    
30.     while t <= 6: # 迭代边界条件
    
31.         y = runge_kutta(y, t, dt, func)
    
32.         t += dt
    
33.         ys.append(y)
    
34.         ts.append(t)
    
35. 
    
36.     plt.plot(ts, ys, label='runge_kutta')
    
37.     plt.legend()
    
38.     plt.show()

复制代码

最佳答案

月排行榜 / 总排行榜

WylLy

2020-6-29 17:42:23

matlab的图中一个x值对应两个y值，但是你写的函数runge_kutta返回值只有一个，所以一个x对应的是一个y值，肯定就跟matlab图不一样了

跳转到最佳答案楼层

WylLy

matlab的图中一个x值对应两个y值，但是你写的函数runge_kutta返回值只有一个，所以一个x对应的是一个y值，肯定就跟matlab图不一样了

kevinheros

没有研究这个的大佬吗？