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

River Chandler

已于 2023-08-28 19:21:08 修改

阅读量1.5k

点赞数 3

分类专栏： # 科学计算原理文章标签：算法数学建模抽象代数 numpy python

于 2023-08-26 19:39:18 首次发布

本文链接：https://blog.csdn.net/Chandler_river/article/details/132515112

版权

科学计算原理专栏收录该内容

24 篇文章 2 订阅

订阅专栏

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

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 方法

$\begin{matrix} y_{n+1}=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)\\ k_1=f(x_n,y_n)\\ k_2=f(x_n+\frac{h}{2},y_n+\frac{hk_1}{2})\\ k_3=f(x_n+\frac{h}{2},y_n+\frac{hk_2}{2})\\ k_4=f(x_n+h,y_i+hk_3) \end{matrix}$

二级四阶隐式 Runge - Kutta 方法

$\begin{matrix} y_{n+1}=y_n+\frac{1}{2}(k_1+k_2)\\ k_1 = f(x_n+\frac{3-\sqrt{3}}{6}h,y_i+\frac{k_1h}{4}+\frac{3-2\sqrt{3}}{12}k_2h)\\ k_2 = f(x_n+\frac{3+\sqrt{3}}{6}h,y_i+\frac{k_2h}{4}+\frac{3+2\sqrt{3}}{12}k_1h) \end{matrix}$

三级六阶隐式 Runge - Kutta 方法

$\begin{matrix} y_{n+1}=y_n+\frac{5}{18}k_1+\frac{4}{9}k_2+\frac{5}{18}k_3\\ k_1=hf(x_n+\frac{5-\sqrt{15}}{10}h,y_n+\frac{5}{36}k_1+\frac{10-3\sqrt{15}}{45}k_2+\frac{25-6\sqrt{15}}{180}k_3)\\ k_2=hf(x_n+\frac{h}{2},y_n+\frac{10+3\sqrt{15}}{72}k_1+\frac{2}{9}k_2+\frac{10-3\sqrt{15}}{12}k_3)\\ k_3=hf(x_n+\frac{5+\sqrt{15}}{10}h,y_n+\frac{25+6\sqrt{15}}{180}k_1+\frac{10+3\sqrt{15}}{45}k_2+\frac{5}{36}k_3) \end{matrix}$

四级四阶 Kutta 方法

$\left\{\begin{matrix} y_{n+1}=y_n + \frac{1}{8}(k_1+8k_2+3k_3+k_4)\\ k_1=hf(x_n,y_n)\\ k_2=hf(x_n+\frac{h}{3},y_n+\frac{k_1}{3})\\ k_3=hf(x_n+\frac{2h}{3},y_n+\frac{k_1}{3}+k_2)\\ k_4=hf(x_n+h,y_n+k_1-k_2+k_3) \end{matrix}\right.$

四级四阶 Milne 方法

$\left\{\begin{matrix} y_{n+1}=y_n + \frac{1}{6}(k_1+(2-\sqrt{2})k_2+(2+\sqrt{2})k_3+k_4)\\ k_1=hf(x_n,y_n)\\ k_2=hf(x_n+\frac{h}{2},y_n+\frac{k_1}{2})\\ k_3=hf(x_n+\frac{h}{2},y_n+\frac{\sqrt{2}-1}{2}k_1+(1-\frac{\sqrt{2}}{2})k_2)\\ k_4=hf(x_n+h,y_n-\frac{\sqrt{2}}{2}k_2+(1+\frac{\sqrt{2}}{2})k_3) \end{matrix}\right.$

求根方法

这个我以后再讲，前面陆陆续续讲过几次，在特征值里面，请大家自行翻阅

常微分方程组的隐式与显式解法

这个就比较酷炫了，因为写起来代码会稍微有一点问题

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)

数值实验

数值实验一

$\begin{pmatrix} y_1'\\ y_2' \end{pmatrix}=\begin{matrix} x^2+y_2^3\\ -y_1^2 \end{matrix}$ $Y_0=\begin{pmatrix} 0\\ 1 \end{pmatrix}$

这一点说明 Milne方法还是有一点不精确的

数值实验二

$\begin{pmatrix} y_1'\\ y_2' \\y_3'\\ y_4' \end{pmatrix}=\begin{matrix} x^2+y_2^3\\ -y_1^2\\y_2-y_4\\y_2y_3 \end{matrix}$ $Y_0=\begin{pmatrix} 0\\ 1\\5 \\7 \end{pmatrix}$

这一点说明 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 秒

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