数值方法2：龙格库塔法在解微分方程中的应用

初学练习，看b站课程，教学为matlab代码，自己写的Python代码，后面会放b站课程链接，感兴趣的同学可以学习观看。

说明：Python初学者，代码可能不够漂亮，欢迎大家批评指正。本系列代码用notebook写的，有些图片显示的命令在pycharm中可能会报错，注释掉就行。

理论部分是课程笔记，有不明白的地方可以看b站课程（一位宝藏老师）。
数值方法2:误差分析，龙格库塔法_哔哩哔哩_bilibilihttps://www.bilibili.com/video/BV1AE411s7e6?spm_id_from=333.999.0.0数值方法2:龙格库塔法，混沌现象_哔哩哔哩_bilibilihttps://www.bilibili.com/video/BV1iE411g7xn?spm_id_from=333.999.0.0

指数函数，这里按照龙格库塔法推导后，又将k1、k2代入、化简，公式变得和二阶欧拉法一样

#!/usr/bin/env python
# coding: utf-8# 数值方法2：龙格库塔法（Runge-Kutta）import numpy as np
import matplotlib.pyplot as plt# 参数
T = 5
N = 100
dt = T/N
t = np.linspace(0,T,N+1)
y0 = 1# 初始化
y = [] # 欧拉一阶
y.append(y0)
# print(y[0])
y1 = []
y1.append(y0) # 欧拉二阶==Runge-Kutta# 求解
for i in range(len(t)-1):y.append(y[i]+y[i]*dt)y1.append(y1[i]+y1[i]*dt+y1[i]*dt*dt/2)# 绘图
plt.plot(t,y,label='Euler')
plt.plot(t,y1,label='Runge-Kutta')# 图例，位置在左上角
plt.plot(t,np.exp(t),label='Real')
plt.legend(loc = 'upper left')plt.title("Runge-Kutta")plt.savefig("Runge-Kutta.png")  # 保存图片plt.show()

指数函数，这里按照四阶龙格库塔法推导后，又将k1~k4代入、化简，公式和欧拉法一样

#!/usr/bin/env python
# coding: utf-8# 数值方法2：龙格库塔法（Runge-Kutta）import numpy as np
import matplotlib.pyplot as plt# 参数
T = 5
N = 100
dt = T/N
t = np.linspace(0,T,N+1)
y0 = 8# 初始化
y = [] # Runge-Kutta2
y.append(y0)y1 = []
y1.append(y0) # Runge-Kutta4# 求解
for i in range(len(t)-1):y.append(y[i]+y[i]*dt+y[i]*dt*dt/2)y1.append(y1[i]+y1[i]*dt+y1[i]*dt*dt/2+y1[i]*dt*dt*dt/6+y1[i]*dt*dt*dt*dt/24)# 绘图
plt.plot(t,y,label='Runge-Kutta2')
plt.plot(t,y1,label='Runge-Kutta4')# 图例，位置在左上角
plt.plot(t,8*np.exp(t),label='Real')
plt.legend(loc = 'upper left')plt.title("Runge-Kutta2vs4")plt.savefig("Runge-Kutta4.png")  # 保存图片plt.show()

洛伦兹系统

#!/usr/bin/env python
# coding: utf-8# 数值方法2：龙格库塔法（Runge-Kutta）
# 洛伦兹系统
# 初值发生一点点改变，随着时间，值会发生很大变化，即对初值敏感，称为混沌现象import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d# 参数
T = 30
N = 10000
dt = T/N
t = []
t = np.linspace(0,T,N+1)
# print(t,len(t))s = 10
b = 8/3
r = 28# 初始条件
y = []
y.append([-1.6,-0.5,21])
# print(y)for n in range(N):k1=[]k2=[]k3=[]k4=[]k1.append(s*(y[n][1]-y[n][0]))k1.append(r*y[n][0]-y[n][1]-y[n][0]*y[n][2])k1.append(-b*y[n][2]+y[n][0]*y[n][1])
#     print(k1)#     y[n][0] => (y[n][0]+k1[0]*dt/2);
#     y[n][1] => (y[n][1]+k1[1]*dt/2);
#     y[n][2] => (y[n][2]+k1[2]*dt/2); k2.append(s*((y[n][1]+k1[1]*dt/2)-y[n][0]+k1[0]*dt/2))k2.append(r*(y[n][0]+k1[0]*dt/2)-(y[n][1]+k1[1]*dt/2)-(y[n][0]+k1[0]*dt/2)*(y[n][2]+k1[2]*dt/2))k2.append(-b*(y[n][2]+k1[2]*dt/2)+(y[n][0]+k1[0]*dt/2)*(y[n][1]+k1[1]*dt/2))
#     print(k2)#     k1 => k2; k3.append(s*((y[n][1]+k2[1]*dt/2)-y[n][0]+k2[0]*dt/2))k3.append(r*(y[n][0]+k2[0]*dt/2)-(y[n][1]+k2[1]*dt/2)-(y[n][0]+k2[0]*dt/2)*(y[n][2]+k2[2]*dt/2))k3.append(-b*(y[n][2]+k2[2]*dt/2)+(y[n][0]+k2[0]*dt/2)*(y[n][1]+k2[1]*dt/2))
#     print(k3)#     k2 => k3; dt/2 => dt;k4.append(s*((y[n][1]+k3[1]*dt)-y[n][0]+k3[0]*dt))k4.append(r*(y[n][0]+k3[0]*dt)-(y[n][1]+k3[1]*dt)-(y[n][0]+k3[0]*dt)*(y[n][2]+k3[2]*dt))k4.append(-b*(y[n][2]+k3[2]*dt)+(y[n][0]+k3[0]*dt)*(y[n][1]+k3[1]*dt))
#     print(k4)y_n = []for m in range(3):y_n.append(y[n][m]+1/6*(k1[m]+2*k2[m]+2*k3[m]+k4[m])*dt)y.append(y_n)
# print(y)Y = np.array(y)
X = np.array(t)
# 绘图
plt.plot(X,Y[:,0])
plt.plot(X,Y[:,1])
plt.plot(X,Y[:,2])plt.title("Lorentz System")plt.savefig("LorentzSystem.png")  # 保存图片
plt.show()#三维
ax = plt.axes(projection='3d')
ax.plot3D(Y[:,0],Y[:,1],Y[:,2])plt.savefig("LorentzSystemPhaseSpace.png")  # 保存图片
plt.show()

数值方法2：龙格库塔法在解微分方程中的应用相关推荐

matlab中函数或变量无法识别怎么办_用MATLAB巧解微分方程实例分析
点"考研竞赛数学"↑可每天"涨姿势"哦! MATLAB巧解微分方程实例分析王少华西安电子科技大学微分方程求解难, 字母一堆看着烦. 写错数字一时爽, 一直 ...
python解常微分方程龙格库_excel实现四阶龙格库塔法runge-kutta解二阶常微分方程范例.xls...
excel实现四阶龙格库塔法runge-kutta解二阶常微分方程范例,rungekutta,四阶rungekutta法,rungekuttamatlab,四阶rungekutta,rungekutt ...
python求解析解,Python解微分方程
Python解微分方程微分方程回顾微分方程:python 解析解(SymPy) 微分方程:python数值解(SciPY) 微分方程组:python数值解微分方程回顾微分方程是用来描述某一类函 ...
线性代数学习笔记7-3：特征值的应用——解微分方程、矩阵的指数函数
之前介绍了求解一阶差分方程,本文介绍求解一阶导数常系数微分方程常系数微分方程的解是指数形式的 e λ t e^{\lambda t} eλt,基于这个事实,我们得以将问题转为线性代数的问题,求解其指 ...
C#，数值计算，解微分方程的龙格-库塔二阶方法与源代码
微分方程含有导数或微分的方程称为微分方程,未知函数为一元函数的微分方程称为常微分方程. 微分方程的阶数微分方程中导数或微分的最高阶数称为微分方程的阶数. 微分方程的解使得微分方程成立的函数称为微 ...
微分方程中的自洽系统（Autonomous system）
微分方程中的自洽系统(Autonomous System) 微分方程中,自洽系统(Autonomous System)表示隐含独立变量的常微分方程系统.特别地,当独立变量是时间 ttt 时,这时的自洽 ...
用python的scipy中的odeint来解常微分方程中的一些细节问题（适用于小白）
用python的scipy中的odeint来解常微分方程中的一些细节问题(适用于小白) 写在前面最近有些需要解决常微分方程的问题,网上查了很多教程都不是很明晰,便自己研究了一段时间,写一点小白初次接 ...
matlab分布傅里叶解微分方程,专题四微分方程的matlab求解
[实验目的及要求] 熟练掌握dsolve,ode45,ode15s语句的调用格式:正确区分符号运算与数值运算的不同,避免出现混淆调用的情况:对实际问题,能够建立微分方程数学模型,能够用dsolve,o ...
详解OpenCV中的Lucas Kanade稀疏光流单应追踪器
详解OpenCV中的Lucas Kanade稀疏光流单应追踪器 1. 效果图 2. 源码参考这篇博客将详细介绍OpenCV中的Lucas Kanade稀疏光流单应追踪器. 光流是由物体或相机的运动 ...

数值方法2：龙格库塔法在解微分方程中的应用

数值方法2：龙格库塔法在解微分方程中的应用相关推荐

最新文章

热门文章