数值方法1：欧拉法在解微分方程中的应用

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

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

理论部分是课程笔记，有不明白的地方可以看b站课程（一位宝藏老师）。

数值方法1:简单的的欧拉法在解微分方程时有多有用（上）_哔哩哔哩_bilibilihttps://www.bilibili.com/video/BV1x7411M7H7?spm_id_from=333.999.0.0

数值方法1:简单的的欧拉法在解微分方程时有多有用（下）:吸引子，极限环，奇怪吸引子_哔哩哔哩_bilibilihttps://www.bilibili.com/video/BV1RE411E7zp?spm_id_from=333.999.0.0

​

 指数函数

#!/usr/bin/env python
# coding: utf-8

# 欧拉法在微分方程中的应用——例一 指数函数

import numpy as np
import matplotlib.pyplot as plt

# 参数
T = 2
N = 100
dt = T/N
t = np.linspace(0,T,N+1)
y0 = 1

# 初始化
y = []
y.append(y0)
# print(y[0])

# 求解
for i in range(len(t)-1):
    y.append(y[i]+y[i]*dt)


# 绘图
plt.plot(t,y,label='Euler')

# 图例，位置在左上角
plt.plot(t,np.exp(t),label='Real')
plt.legend(loc = 'upper left')

plt.title("ExponentialFunction")

plt.savefig("ExponentialFunction.png")  # 保存图片

plt.show()

 ​

谐振子 

#!/usr/bin/env python
# coding: utf-8

# 欧拉法解二阶微分方程——例二 谐振子

import numpy as np
import matplotlib.pyplot as plt

# 参数
T = 20
N = 10000
dt = T/N
t = np.linspace(0,T,N+1)
k = 1
m = 1
y10 = 1
y20 = 1

# 初始化
y1 = []
y1.append(y10)
# print(y[0])
y2 = []
y2.append(y20)
# print(y[0])

# 求解
for i in range(len(t)-1):
    y1.append(y1[i]+y2[i]*dt)
    y2.append(y2[i]-k/m*y1[i]*dt)

# 绘图
plt.plot(t,y1,label='y1')
plt.plot(t,y2,label='y2')

# 图例，位置在左上角
plt.legend(loc = 'upper left')

plt.title("Harmonic Oscillator")

plt.savefig("HarmonicOscillator.png")  # 保存图片

plt.show()

# 绘图,相空间
plt.plot(y1,y2)
plt.title("Phase Space")

plt.savefig("PhaseSpace.png")  # 保存图片

plt.show()

​ ​

阻尼谐振子 

#!/usr/bin/env python
# coding: utf-8

# 欧拉法解二阶微分方程——例3 阻尼谐振子
# 形成吸引子（attractor）

import numpy as np
import matplotlib.pyplot as plt

# 参数
T = 20
N = 10000
dt = T/N
t = np.linspace(0,T,N+1)
k = 1
m = 1
g = 0.5 # 阻尼
y10 = 1
y20 = 1

# 初始化
y1 = []
y1.append(y10)
# print(y[0])
y2 = []
y2.append(y20)
# print(y[0])

# 求解
for i in range(len(t)-1):
    y1.append(y1[i]+y2[i]*dt)
    y2.append(y2[i]-(k/m*y1[i]+g*y2[i])*dt)

# 绘图
plt.plot(t,y1,label='y1')
plt.plot(t,y2,label='y2')

# 图例，位置在左上角
plt.legend(loc = 'upper left')
plt.title("Damped Oscillators")
plt.savefig("DampedOscillators.png")  # 保存图片
plt.show()

# 绘图,相空间
plt.plot(y1,y2)
plt.title("Attractor")

plt.savefig("Attractor.png")  # 保存图片

plt.show()

​​

单摆 

#!/usr/bin/env python
# coding: utf-8

# 欧拉法解二阶微分方程——例4 单摆

import numpy as np
import matplotlib.pyplot as plt

# 参数
T = 6
N = 10000
dt = T/N
t = np.linspace(0,T,N+1)
g = 9.8
l = 1

# 初始化
y1 = []
y1.append(3) # 初始位移
# y1.append(np.pi) # 初始位移
y2 = []
y2.append(0) # 初始速度
# y2.append(0.01) # 初始速度

# 求解
for i in range(len(t)-1):
    y1.append(y1[i]+y2[i]*dt)
    y2.append(y2[i]-g/l*np.sin(y1[i])*dt)

# 绘图
plt.plot(t,y1,label='y1')
plt.plot(t,y2,label='y2')

# 图例，位置在左上角
plt.legend(loc = 'upper left')
plt.title("simple pendulum")
plt.savefig("SimplePendulum1.png")  # 保存图片
plt.show()

# 绘图,相空间
plt.plot(y1,y2)
plt.title("Phase Space")
# plt.xlim(-3.0,3.0)
# plt.xlim(-3.0,3.0)
plt.savefig("SimplePendulumPhaseSpace1.png")  # 保存图片

plt.show()

 ​​

 瑞利方程——极限环

#!/usr/bin/env python
# coding: utf-8

# 欧拉法解二阶微分方程——例5 瑞利方程（Rayleigh equation）
# 形成极限环

import numpy as np
import matplotlib.pyplot as plt

# 参数
T = 20
N = 10000
dt = T/N
t = np.linspace(0,T,N+1)
e = 0.1

# 初始化
y1 = []
y1.append(5) # 初始位移

y2 = []
y2.append(0) # 初始速度

# 求解
for i in range(len(t)-1):
    y1.append(y1[i]+y2[i]*dt)
    y2.append(y2[i]+(e*y2[i]*(1-np.power(y2[i],2))-y1[i])*dt)

# 绘图
plt.plot(t,y1,label='y1')
plt.plot(t,y2,label='y2')

# 图例，位置在左上角
plt.legend(loc = 'upper left')
plt.title("Rayleigh equation")
plt.savefig("Rayleighequation.png")  # 保存图片
plt.show()

# 绘图,相空间
plt.plot(y1,y2)
plt.title("Limit Cycle")
# plt.xlim(-3.0,3.0)
# plt.xlim(-3.0,3.0)
plt.savefig("LimitCycle.png")  # 保存图片

plt.show()

 ​​

 洛伦兹系统——奇怪吸引子

#!/usr/bin/env python
# coding: utf-8

# 欧拉法解二阶微分方程——例6 洛伦系统
# 形成奇怪吸引子（StrangeAttractor）

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 参数
T = 50
N = 10000
dt = T/N
t = np.linspace(0,T,N+1)

sigam = 10
beita = 8/3
rou = 21

# 初始化
y1 = []
y1.append(1) # 初始位移

y2 = []
y2.append(0) # 初始速度

y3 = []
y3.append(0) # 初始速度


# 求解
for i in range(len(t)-1):
    y1.append(y1[i]+sigam*(y2[i]-y1[i])*dt)
    y2.append(y2[i]+(rou*y1[i]-y2[i]-y1[i]*y3[i])*dt)
    y3.append(y3[i]+(-beita*y3[i]+y1[i]*y2[i])*dt)

# 3D曲面
get_ipython().run_line_magic('matplotlib', 'notebook') # notebook显示图片，pycharm中注释掉
#定义坐标轴
fig = plt.figure()
ax1 = plt.axes(projection='3d')
ax1.plot3D(y1,y2,y3)    #绘制空间曲线
plt.savefig("StrangeAttractor.png")  # 保存图片
plt.show()

 ​