常微分方程的龙格库塔显式与隐式解法

• 好习惯，讲问题之前先来介绍一下最近生活状况。
  • 得到了 熟人的 熟人认证 很好 很荣幸了 属于是
• 先上全部代码与效果图

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolveclass odesolver():def __init__(self, f, X_start=0, X_end=1, Y_start=1, dY_start=1, h=0.01):self.f = fself.h = hself.X = np.arange(X_start, X_end, self.h)self.Y = np.zeros(self.X.size)self.Y[0] = Y_startself.Y[1] = Y_start + self.h * dY_startself.tol = 1e-6def __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.Ydef 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.Ydef 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.Ydef 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.Ydef 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.Yc = 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 fsolveclass 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 = fself.h = hself.X = np.arange(X_start, X_end, self.h)self.n = Y_start.sizeself.Y = np.zeros((self.n, self.X.size))#第一个参数表示元 第二个参数表示变量self.Y[:, 0] = Y_startself.Y[:, 1] = Y_start + self.h * dY_startself.tol = 1e-6def __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.Ydef 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.Ydef 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.Ydef 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)/6f1_y = self.Y[:, i-1]+k1/4*self.h+(3-2*3**0.5)/12*k2*self.hf1_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)/6f2_y = self.Y[:, i-1]+k2/4*self.h+(3+2*3**0.5)/12*k1*self.hf2_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 docsol = 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.Ydef 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)/10f1_y = self.Y[:, i-1]+k1*5/36+(10-3*15**0.5)/45*k2+(25-6*15**0.5)/180*k3f1_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/2f2_y = self.Y[:, i-1]+k1*(10+3*15**0.5)/72+2/9*k2+(10-3*15**0.5)/12*k3f2_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)/10f3_y = self.Y[:, i-1]+k1*(25+6*15**0.5)/180+(10+3*15**0.5)/45*k2+5/36*k3f3_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 docsol = 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.Ydef 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)/10f1_y = self.Y[:, i-1]+k1*5/36+(10-3*15**0.5)/45*k2+(25-6*15**0.5)/180*k3f1_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/2f2_y = self.Y[:, i-1]+k1*(10+3*15**0.5)/72+2/9*k2+(10-3*15**0.5)/12*k3f2_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)/10f3_y = self.Y[:, i-1]+k1*(25+6*15**0.5)/180+(10+3*15**0.5)/45*k2+5/36*k3f3_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 docsol = 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 = paramreturn 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 datetimec = 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 秒

• 记得给电脑降温，这很重要

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！