鱼C论坛

 找回密码
 立即注册
查看: 152|回复: 2

微分方程的龙格-库塔法Python编程请教

[复制链接]
最佳答案
1 
发表于 2020-6-29 17:42:22 | 显示全部楼层 |阅读模式
10鱼币
有个微分方程要用龙格-库塔法求解,但我编程求解后的图和实际的图不一样,是否有了解的鱼油看看问题所在?多谢!!!
1. 原微分方程:
01.jpg
2.边界条件:(其中,h=8mm,理解为当x=0时,y=8)
02.jpg
3.别人用Matlab编程得到的曲线图:
03.jpg
4.我用Python编写的程序及对应的曲线图:
  1. import numpy as np
  2. import matplotlib.pyplot as plt

  3. # 龙格-库塔法函数定义
  4. def runge_kutta(y, x, dx, f):
  5.     """ y is the initial value for y
  6.         x is the initial value for x
  7.         dx is the time step in x
  8.         f is derivative of function y(t)
  9.     """
  10.     k1 = dx * f(y, t)
  11.     k2 = dx * f(y + 0.5 * k1, x + 0.5 * dx)
  12.     k3 = dx * f(y + 0.5 * k2, x + 0.5 * dx)
  13.     k4 = dx * f(y + k3, x + dx)
  14.     return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6.

  15. if __name__=='__main__':
  16.     t = 0.0
  17.     y = 8.0 # 初始值为(0,8)
  18.     dt = 0.000001 #步长
  19.     ys, ts = [], []

  20.     # 微分方程函数定义
  21.     def func(y, t):
  22.         v1 = np.sqrt(t**2 + y**2)
  23.         v2 = 1000 * t / v1
  24.         v3 = np.sqrt((v2 - t)**2 + (1000 - y)**2)
  25.         return (v1*(v2-t)-1.49*t*v3)/(1.49*y*v3-(1000-y)*v1)
  26.    
  27.     while t <= 6: # 迭代边界条件
  28.         y = runge_kutta(y, t, dt, func)
  29.         t += dt
  30.         ys.append(y)
  31.         ts.append(t)

  32.     plt.plot(ts, ys, label='runge_kutta')
  33.     plt.legend()
  34.     plt.show()
复制代码

Figure_1.png

想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
最佳答案
1 
 楼主| 发表于 7 天前 | 显示全部楼层
没有研究这个的大佬吗?
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
最佳答案
8 
发表于 7 天前 | 显示全部楼层
matlab的图中一个x值对应两个y值,但是你写的函数runge_kutta返回值只有一个,所以一个x对应的是一个y值,肯定就跟matlab图不一样了
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

关闭

小甲鱼强烈推荐上一条 /1 下一条

小黑屋|手机版|Archiver|鱼C工作室 ( 粤ICP备18085999号-1

GMT+8, 2020-7-9 17:11

Powered by Discuz! X3.4

Copyright © 2001-2020, Tencent Cloud.

快速回复 返回顶部 返回列表