数值方法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()