研究生课程《应用数值分析》结课了，使用python简单记录了下常微分方程数值解法。

• 2022.11.26 Update: 文末补充C语言实现（C11标准）

向前欧拉法

$$\{\begin{array}{l}{y}_{i+1}=y_i+h_{i}f(x_i,y_i)\\ y_0=y(a)\end{array}$$

from pylab import *
import warnings
warnings.filterwarnings('ignore')

求解初值问题
$$\{\begin{array}{lr}{y}^{\prime }=x-y+1& 0\le x\le 1\\ y(0)=1& \end{array}$$

该微分方程的精确解为：$$y=x+exp(-x)$$

def f(t,y):
    '''
    求解的微分方程，
    '''
    return t-y+1

def euler_forward(f,a=0,b=1,ya=1,h=0.1,verbose=True):
    '''向前欧拉法
    Args
    ----------
    f： callable function
        需要求解的函数
    a: float
        求解区间起始值
    b：float
        求解区间终止值
    ya：float
        起始条件，ya=y(a)
    h：float
        求解步长（区间[a,b]n等分）
    verbose：logical，default is True
        显示迭代结果
    
    Returns
    ----------
    res：list like
        返回向前欧拉发求解的结果
    '''
#     i = 0
    res = []
    xi = a 
    yi = ya
    while xi<=b: # 在求解区间范围
        y = yi + h*f(xi,yi)
        if verbose:
            print('xi:{:.2f}, yi:{:.6f}'.format(xi,yi))
        res.append(y)
        xi, yi = xi+h, y
    
    return res

res = euler_forward(f,a=0,b=1,ya=1,h=0.1,verbose=True)

xi:0.00, yi:1.000000
xi:0.10, yi:1.000000
xi:0.20, yi:1.010000
xi:0.30, yi:1.029000
xi:0.40, yi:1.056100
xi:0.50, yi:1.090490
xi:0.60, yi:1.131441
xi:0.70, yi:1.178297
xi:0.80, yi:1.230467
xi:0.90, yi:1.287420
xi:1.00, yi:1.348678

龙格库塔(Runge-Kutta methods)

三阶龙格库塔法

$$\{\begin{array}{lr}{y}_{i+1}=y_i+\frac{1}{6}(k_1+2k_2+k_3)& \\ k_1=hf(x_i,y_i)& \\ k_2=hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}{k}_1)& \\ k_3=hf(x_i+h,y_i-k_1+2k_2)& \\ y_0=y(a),i=0,1,\cdots ,n-1& \end{array}$$

四阶龙格库塔法

$$\{\begin{array}{lr}{y}_{i+1}=y_i+\frac{1}{6}(k_1+2k_2+2k_3+k_4)& \\ k_1=hf(x_i,y_i)& \\ k_2=hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}{k}_1)& \\ k_3=hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}{k}_2)& \\ k_4=hf(x_i+h,y_i+k_3)& \\ y_0=y(a),i=0,1,\cdots ,n-1& \end{array}$$

def runge_kutta(f,a=0,b=1,ya=1,h=0.1,verbose=True):
    '''四阶龙格库塔法
    Args
    ----------
    f： callable function
        需要求解的函数
    a: float
        求解区间起始值
    b：float
        求解区间终止值
    ya：float
        起始条件，ya=y(a)
    h：float
        求解步长（区间[a,b]n等分）
    verbose：logical，default is True
        显示迭代结果
    
    Returns
    ----------
    res：list like
        返回向前欧拉发求解的结果
    '''
    res = []
    xi = a 
    yi = ya
    
    while xi <= b: # 在求解区间范围
        k1 = h * f(xi, yi)
        k2 = h * f(xi + h/2, yi + k1/2)
        k3 = h * f(xi + h/2, yi + k2/2)
        k4 = h * f(xi + h, yi + k3)
        y = yi + 1/6 * (k1 + 2*k2 + 2*k3 + k4)
        if verbose:
            print('xi:{:.2f}, yi:{:.10f}'.format(xi,yi))
        res.append(y)
        xi, yi = xi+h, y
    
    return res
res = runge_kutta(f,a=0,b=1,ya=1,h=0.1,verbose=True)

xi:0.00, yi:1.0000000000
xi:0.10, yi:1.0048375000
xi:0.20, yi:1.0187309014
xi:0.30, yi:1.0408184220
xi:0.40, yi:1.0703202889
xi:0.50, yi:1.1065309344
xi:0.60, yi:1.1488119344
xi:0.70, yi:1.1965856187
xi:0.80, yi:1.2493292897
xi:0.90, yi:1.3065699912
xi:1.00, yi:1.3678797744

改进欧拉法

$$\{\begin{array}{lr}{y}_p=y_i+hf(x_i,y_i)& \\ y_{i+1}=y_i+\frac{h}{2}[f(x_i,y_i)+f(x_i,y_p)]& i=0,1,\cdots ,n-1\\ y_0=y(a)& \end{array}$$

求解初值问题
$$\{\begin{array}{lr}{y}^{\prime }=-y& (0\le x\le 1)\\ y(0)=1& \end{array}$$
精确解为$$y=exp(-x)$$

def f(t,y):
    '''
    精确解为y=exp(-x)
    '''
    return -y

def improved_euler(f,a=0,b=1,ya=1,h=0.1,verbose=True):
    '''改进欧拉法
    Args
    ----------
    f： callable function
        需要求解的函数
    a: float
        求解区间起始值
    b：float
        求解区间终止值
    ya：float
        起始条件，ya=y(a)
    h：float
        求解步长（区间[a,b]n等分）
    verbose：logical，default is True
        显示迭代结果
    
    Returns
    ----------
    res：list like
        返回向前欧拉发求解的结果
    '''
    res = []
    xi = a 
    yi = ya
    
    while xi <= b: # 在求解区间范围
        yp = yi + h*f(xi, yi)
        y = yi + h/2 * (f(xi, yi) + f(xi, yp))
        if verbose:
            print('xi:{:.2f}, yi:{:.6f}'.format(xi,yi))
        res.append(y)
        xi, yi = xi+h, y
    
    return res
res = improved_euler(f,a=0,b=1,ya=1,h=0.1,verbose=True)

xi:0.00, yi:1.000000
xi:0.10, yi:0.905000
xi:0.20, yi:0.819025
xi:0.30, yi:0.741218
xi:0.40, yi:0.670802
xi:0.50, yi:0.607076
xi:0.60, yi:0.549404
xi:0.70, yi:0.497210
xi:0.80, yi:0.449975
xi:0.90, yi:0.407228
xi:1.00, yi:0.368541

C语言实现（C11标准）:

/* C11 standard*/
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <stdbool.h>

// ode 1: dy/dx = x - y + 1; y(0)=1
// exac solution: y = x + exp(-x)
float eq1(float x, float y)
{
    return x - y + 1.0;
}

float f(float x)
{
    return x + exp(-x);
}

void euler_forward(float (*func)(float, float), float x0, float y0, float h, float x, bool verb)
{
    int iter = 0;
    while (x0 < x)
    {
        iter++;
        y0 = y0 + h * func(x0, y0);
        x0 = x0 + h;
        if (verb) {
            printf("Iter: %d, xi: %f, yi: %f, yi(true): %f\n", iter, x0, y0, f(x0));
        }
    }
    printf("-----------------------------------------------\n");
}

int runge_kutta(float (*func)(float, float), float x1, float x2, float y0, float h, bool verb)
{
    assert (x1 < x2);
    float k1, k2, k3, k4;
    float y = y0;
    int iter = 0;
    while (x1 < x2)
    {
        iter ++ ;
        k1 = h * func(x1, y0);
        k2 = h * func(x1+h/2., y0+k1/2);
        k3 = h * func(x1+h/2., y0+k2/2);
        k4 = h * func(x1 + h, y0 + k3);
        y = y0 + 1. / 6. * (k1 + 2 * k2 + 2 * k3 + k4);
        x1 = x1 + h;
        y0 = y;
        if (verb) 
        {
            // printf("k1 = %f, k2 = %f, k3 = %f, k4 = %f\n", k1, k2, k3, k4);
            printf("Iter: %d, xi: %f, yi: %f, yi(true): %f\n", iter, x1, y, f(x1));
        }
    }
    printf("--------------------------------------------\n");
    return 0;
}
int main()
{
    float x0 = 0; // initial x
    float y0 = 1; // initial y
    float h = 0.1;
    float x = 1;
    // [1] euler_forward method
    euler_forward(eq1, x0, y0, h, x, true);

    // [2] RK4 method
    runge_kutta(eq1, x0, x, y0, h, true);
    return 0;
}