- 好习惯,讲问题之前先来介绍一下最近生活状况。
- 得到了 熟人的 熟人认证 很好 很荣幸了 属于是
- 先上全部代码与效果图
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolve
class odesolver():
def __init__(self, f, X_start=0, X_end=1, Y_start=1, dY_start=1, h=0.01):
self.f = f
self.h = h
self.X = np.arange(X_start, X_end, self.h)
self.Y = np.zeros(self.X.size)
self.Y[0] = Y_start
self.Y[1] = Y_start + self.h * dY_start
self.tol = 1e-6
def __str__(self):
return f"y'(x) = f(x) = ({self.f}) only one variable"
def RK4(self):
for i in range(1, self.X.size):
k1 = self.f(self.X[i-1], self.Y[i-1])
k2 = self.f(self.X[i-1]+self.h/2, self.Y[i-1]+1/2*self.h*k1)
k3 = self.f(self.X[i-1]+self.h/2, self.Y[i-1]+1/2*self.h*k2)
k4 = self.f(self.X[i-1]+self.h, self.Y[i-1]+self.h*k3)
self.Y[i] = self.Y[i-1] + 1/6 * self.h * (k1 + 2*k2 + 2*k3 + k4)
return self.Y
def IRK4(self):
for i in range(1, self.X.size):
f1 = lambda k:self.f(self.X[i-1]+self.h*(3-3**0.5)/6, \
self.Y[i-1]+k[0]/4*self.h+(3-2*3**0.5)/12*k[1]*self.h) - k[0]
f2 = lambda k:self.f(self.X[i-1]+self.h*(3+3**0.5)/6, \
self.Y[i-1]+k[1]/4*self.h+(3+2*3**0.5)/12*k[0]*self.h) - k[1]
sol = fsolve(lambda k: [f1(k), f2(k)], [0, 0])
self.Y[i] = self.Y[i-1] + 1/2 * self.h * (sol[0] + sol[1])
return self.Y
def IRK6(self):
for i in range(1, self.X.size):
f1 = lambda k:self.h * self.f(self.X[i-1]+self.h*(5-15**0.5)/10, \
self.Y[i-1]+k[0]*5/36 +(10-3*15**0.5)/45*k[1] +(25-6*15**0.5)/180*k[2]) - k[0]
f2 = lambda k:self.h * self.f(self.X[i-1]+self.h/2, \
self.Y[i-1]+k[0]*(10+3*15**0.5)/72 +2/9*k[1] +(10-3*15**0.5)/12*k[2]) - k[1]
f3 = lambda k:self.h * self.f(self.X[i-1]+self.h*(5+15**0.5)/10, \
self.Y[i-1]+k[0]*(25+6*15**0.5)/180 +(10+3*15**0.5)/45*k[1] +5/36*k[2]) - k[2]
sol = fsolve(lambda k: [f1(k), f2(k), f3(k)], [0, 0, 0])
self.Y[i] = self.Y[i-1] + (5/18*sol[0] + 4/9*sol[1] + 5/18*sol[2])
return self.Y
def Milne(self):
for i in range(1, self.X.size):
k1 = self.h * self.f(self.X[i-1], self.Y[i-1])
k2 = self.h * self.f(self.X[i-1]+self.h/3, self.Y[i-1]+1/3*k1)
k3 = self.h * self.f(self.X[i-1]+self.h*2/3, self.Y[i-1]+1/3*k1+k2)
k4 = self.h * self.f(self.X[i-1]+self.h, self.Y[i-1]+k1-k2+k3)
self.Y[i] = self.Y[i-1] + 1/8 * (k1 + 8*k2 + 3*k3 + k4)
return self.Y
def Kutta(self):
for i in range(1, self.X.size):
k1 = self.h * self.f(self.X[i-1], self.Y[i-1])
k2 = self.h * self.f(self.X[i-1]+self.h/2, self.Y[i-1]+1/2*k1)
k3 = self.h * self.f(self.X[i-1]+self.h/2, self.Y[i-1]+(2**0.5-1)/2*k1+(1-2**0.5/2)*k2)
k4 = self.h * self.f(self.X[i-1]+self.h, self.Y[i-1]-2**0.5/2*k2+(1+2**0.5/2)*k3)
self.Y[i] = self.Y[i-1] + 1/6 * (k1 + (2-2**0.5)*k2 + (2+2**0.5)*k3 + k4)
return self.Y
c = odesolver(f=lambda x,y:x*y)
x = np.arange(0, 1, 0.01)
y1 = c.RK4()
y2 = c.IRK4()
y3 = c.IRK6()
y4 = c.Milne()
y5 = c.Kutta()
f_true = lambda x:np.exp(x**2/2)
y_true = f_true(x)
plt.plot(x, y1, label="RK4")
plt.plot(x, y2, label="IRK4")
plt.plot(x, y3, label="IRK6")
plt.plot(x, y4, label="Milne")
plt.plot(x, y5, label="Kutta")
plt.plot(x, y_true, label="true")
plt.legend()
plt.pause(0.01)
- 你就说这个代码写得规不规范,整不整齐
- 效果图
常微分方程的隐式与显式解法详解
- 说到这一点,就必须将一个故事
- 有天 博主在网上逛 但是发现这里竟然没有一篇能完整得给出解微分方程隐式方法的完整代码的博客,于是决定写一篇
- 什么是显式方法呢,就是每一步的值仅有上一步或前几步有关
- 什么是隐式方法呢,就是对于每一步必须要求解一个方程以得到该步的近似值
- 很好,我觉得我解释得极其精确无误了
常用公式
四级四阶显式 Runge - Kutta 方法
二级四阶隐式 Runge - Kutta 方法
三级六阶隐式 Runge - Kutta 方法
四级四阶 Kutta 方法
四级四阶 Milne 方法
求根方法
- 这个我以后再讲,前面陆陆续续讲过几次,在特征值里面,请大家自行翻阅
常微分方程组的隐式与显式解法
- 这个就比较酷炫了,因为写起来代码会稍微有一点问题
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolve
class odessolver():
def __init__(self, f, Y_start=np.array([0, 1]), dY_start=np.array([0, 0]), \
X_start=0, X_end=1, h=0.01):
self.f = f
self.h = h
self.X = np.arange(X_start, X_end, self.h)
self.n = Y_start.size
self.Y = np.zeros((self.n, self.X.size))
#第一个参数表示元 第二个参数表示变量
self.Y[:, 0] = Y_start
self.Y[:, 1] = Y_start + self.h * dY_start
self.tol = 1e-6
def __str__(self):
return f"y'(x) = f(x) = ({self.f}) variables"
def RK4(self):
for i in range(1, self.X.size):
k1 = self.f(self.X[i-1] , self.Y[:, i-1])
k2 = self.f(self.X[i-1] +self.h/2 , self.Y[:, i-1]+1/2*self.h*k1)
k3 = self.f(self.X[i-1] +self.h/2 , self.Y[:, i-1]+1/2*self.h*k2)
k4 = self.f(self.X[i-1] +self.h , self.Y[:, i-1]+ self.h*k3)
self.Y[:, i] = self.Y[:, i-1] +self.h/6 * (k1 + 2*k2 + 2*k3 + k4)
return self.Y
def Milne(self):
for i in range(1, self.X.size):
k1 = self.h * self.f(self.X[i-1] , self.Y[:, i-1])
k2 = self.h * self.f(self.X[i-1] +self.h/3 , self.Y[:, i-1] + 1/3*k1)
k3 = self.h * self.f(self.X[i-1] +self.h*2/3 , self.Y[:, i-1] + 1/3*k1+k2)
k4 = self.h * self.f(self.X[i-1] +self.h , self.Y[:, i-1] + k1 - k2 + k3)
self.Y[:, i] = self.Y[:, i-1] + 1/8 * (k1 + 8*k2 + 3*k3 + k4)
return self.Y
def Kutta(self):
for i in range(1, self.X.size):
k1 = self.h * self.f(self.X[i-1] , self.Y[:, i-1])
k2 = self.h * self.f(self.X[i-1]+self.h/2 , self.Y[:, i-1] + k1/2)
k3 = self.h * self.f(self.X[i-1]+self.h/2 , self.Y[:, i-1] + (2**0.5-1)/2*k1+(1-2**0.5/2)*k2)
k4 = self.h * self.f(self.X[i-1]+self.h , self.Y[:, i-1] - 2**0.5/2*k2+(1+2**0.5/2)*k3)
self.Y[:, i] = self.Y[:, i-1] + 1/6 * (k1 + (2-2**0.5)*k2 + (2+2**0.5)*k3 + k4)
return self.Y
def IRK4(self):
for i in range(1, self.X.size):
def f1(k1, k2):
f1_x = self.X[i-1] + self.h*(3-3**0.5)/6
f1_y = self.Y[:, i-1]+k1/4*self.h+(3-2*3**0.5)/12*k2*self.h
f1_res = self.f(f1_x, f1_y)
return np.array([f1_res[i] for i in range(self.n)])
def f2(k1, k2):
f2_x = self.X[i-1] + self.h*(3+3**0.5)/6
f2_y = self.Y[:, i-1]+k2/4*self.h+(3+2*3**0.5)/12*k1*self.h
f2_res = self.f(f2_x, f2_y)
return np.array([f2_res[i] for i in range(self.n)])
def func(k):
k1 = np.array([k[i] for i in range(self.n)])
k2 = np.array([k[i+self.n] for i in range(self.n)])
doc = []
for i in range(self.n):
doc.append((k1 - f1(k1, k2))[i])
for i in range(self.n):
doc.append((k2 - f2(k1, k2))[i])
return doc
sol = fsolve(func, np.zeros(self.n*2))
self.Y[:, i] = self.Y[:, i-1] + 1/2 * self.h * (sol[:self.n] + sol[self.n:])
return self.Y
def IRK6(self):
for i in range(1, self.X.size):
def f1(k1, k2, k3):
f1_x = self.X[i-1] + self.h*(5-15**0.5)/10
f1_y = self.Y[:, i-1]+k1*5/36+(10-3*15**0.5)/45*k2+(25-6*15**0.5)/180*k3
f1_res = self.f(f1_x, f1_y)
return np.array([f1_res[i] for i in range(self.n)])
def f2(k1, k2, k3):
f2_x = self.X[i-1] + self.h/2
f2_y = self.Y[:, i-1]+k1*(10+3*15**0.5)/72+2/9*k2+(10-3*15**0.5)/12*k3
f2_res = self.f(f2_x, f2_y)
return np.array([f2_res[i] for i in range(self.n)])
def f3(k1, k2, k3):
f3_x = self.X[i-1] + self.h*(5+15**0.5)/10
f3_y = self.Y[:, i-1]+k1*(25+6*15**0.5)/180+(10+3*15**0.5)/45*k2+5/36*k3
f3_res = self.f(f3_x, f3_y)
return np.array([f3_res[i] for i in range(self.n)])
def func(k):
k1 = np.array([k[i] for i in range(self.n)])
k2 = np.array([k[i+self.n] for i in range(self.n)])
k3 = np.array([k[i+2*self.n] for i in range(self.n)])
doc = []
for i in range(self.n):
doc.append((k1 - self.h * f1(k1, k2, k3))[i])
for i in range(self.n):
doc.append((k2 - self.h * f2(k1, k2, k3))[i])
for i in range(self.n):
doc.append((k3 - self.h * f3(k1, k2, k3))[i])
return doc
sol = fsolve(func, np.zeros(self.n*3))
self.Y[:, i] = self.Y[:, i-1] +\
5/18 * sol[:self.n] +\
4/9 * sol[self.n:2*self.n] +\
5/18 * sol[2*self.n:]
return self.Y
def test_fun(x, Y):
return np.array([x**2+Y[1]**3, -Y[0]**2])
c = odessolver(test_fun)
x = np.arange(0, 1, 0.01)
y1 = c.Kutta()
plt.plot(x, y1[0, :], label="Kutta1")
plt.plot(x, y1[1, :], label="Kutta2")
y2 = c.Milne()
plt.plot(x, y1[0, :], label="Milne1")
plt.plot(x, y1[1, :], label="Milne2")
y3 = c.RK4()
x = np.arange(0, 1, 0.01)
plt.plot(x, y1[0, :], label="RK41")
plt.plot(x, y1[1, :], label="RK42")
y4 = c.IRK4()
x = np.arange(0, 1, 0.01)
plt.plot(x, y4[0, :], label="IRK41")
plt.plot(x, y4[1, :], label="IRK42")
y5 = c.IRK6()
x = np.arange(0, 1, 0.01)
plt.plot(x, y5[0, :], label="IRK61")
plt.plot(x, y5[1, :], label="IRK62")
plt.legend()
plt.pause(0.01)
数值实验
数值实验一
- 这一点说明 Milne方法还是有一点不精确的
数值实验二
- 这一点说明 RK4 IRK4 IRK6 效果都十分牛
- Kutta 方法 和 Milne 方法 效果就不好了
- 总体概览图
含参常微分方程组的求解
- 数学模型中我们会常见到带参数的微分方程组,以三级六阶龙格库塔方法为例,我们对这玩意进行一下小小的微操
def IRK6(self, Param=None):
for i in range(1, self.X.size):
def f1(k1, k2, k3):
f1_x = self.X[i-1] + self.h*(5-15**0.5)/10
f1_y = self.Y[:, i-1]+k1*5/36+(10-3*15**0.5)/45*k2+(25-6*15**0.5)/180*k3
f1_res = self.f(f1_x, f1_y, param = Param)
return np.array([f1_res[i] for i in range(self.n)])
def f2(k1, k2, k3):
f2_x = self.X[i-1] + self.h/2
f2_y = self.Y[:, i-1]+k1*(10+3*15**0.5)/72+2/9*k2+(10-3*15**0.5)/12*k3
f2_res = self.f(f2_x, f2_y, param = Param)
return np.array([f2_res[i] for i in range(self.n)])
def f3(k1, k2, k3):
f3_x = self.X[i-1] + self.h*(5+15**0.5)/10
f3_y = self.Y[:, i-1]+k1*(25+6*15**0.5)/180+(10+3*15**0.5)/45*k2+5/36*k3
f3_res = self.f(f3_x, f3_y, param = Param)
return np.array([f3_res[i] for i in range(self.n)])
def func(k):
k1 = np.array([k[i] for i in range(self.n)])
k2 = np.array([k[i+self.n] for i in range(self.n)])
k3 = np.array([k[i+2*self.n] for i in range(self.n)])
doc = []
for i in range(self.n):
doc.append((k1 - self.h * f1(k1, k2, k3))[i])
for i in range(self.n):
doc.append((k2 - self.h * f2(k1, k2, k3))[i])
for i in range(self.n):
doc.append((k3 - self.h * f3(k1, k2, k3))[i])
return doc
sol = fsolve(func, np.zeros(self.n*3))
self.Y[:, i] = self.Y[:, i-1] +\
5/18 * sol[:self.n] +\
4/9 * sol[self.n:2*self.n] +\
5/18 * sol[2*self.n:]
return self.Y
- 而后
def test_fun(x, Y, param=[1, 1, 1]):
a, b, c = param
return np.array([a*x**2+b*Y[1]**3, -c*Y[0]**2])
##c = odessolver(test_fun)
含参常微分方程组求解的并行加速
from paraodessolver import *
import multiprocessing as mp
import datetime
c = odessolver(test_fun)
def final_fun(name,para):
for num in para:
result = c.IRK6(Param=num)
return {name:result}
if __name__ == '__main__':
start_time = datetime.datetime.now()
num_cores = int(mp.cpu_count())
pool = mp.Pool(num_cores)
param_dict = {'task1': [[1, 2, 4],[2, 4, 5],[3, 4, 1]],
'task2': [[3, 3, 3],[2, 2, 2],[1, 1, 6]],
'task3': [[8, 6, 4],[3, 5, 6],[1, 2, 8]],
'task4': [[9, 8, 5],[6, 6, 7],[3, 6, 9]],
'task5': [[1, 2, 4],[2, 4, 5],[3, 4, 1]],
'task6': [[3, 3, 3],[2, 2, 2],[1, 1, 6]],
'task7': [[8, 6, 4],[3, 5, 6],[1, 2, 8]],
'task8': [[9, 8, 5],[6, 6, 7],[3, 6, 9]]}
results = [pool.apply_async(final_fun, args=(name, param)) for name, param in param_dict.items()]
results = [p.get() for p in results]
end_time = datetime.datetime.now()
use_time = (end_time - start_time).total_seconds()
print("多进程计算 共消耗: " + "{:.2f}".format(use_time) + " 秒")
## print(results)
start_time = datetime.datetime.now()
for i in range(3):
c.IRK6(Param=[1, 2, 3])
end_time = datetime.datetime.now()
print("单进程计算 共消耗: " + "{:.2f}".format(use_time) + " 秒")
多进程计算 共消耗: 2.26 秒
单进程计算 共消耗: 2.26 秒
- 效率提高了区区八倍而已,并不是很强的啦
多进程计算 共消耗: 2.16 秒
单进程计算 共消耗: 2.16 秒
多进程计算 共消耗: 2.68 秒
单进程计算 共消耗: 2.68 秒
多进程计算 共消耗: 2.58 秒
单进程计算 共消耗: 2.58 秒
多进程计算 共消耗: 2.55 秒
单进程计算 共消耗: 2.55 秒
多进程计算 共消耗: 2.53 秒
单进程计算 共消耗: 2.53 秒
- 记得给电脑降温,这很重要