研究生课程《应用数值分析》结课了,使用python简单记录了下常微分方程数值解法。
- 2022.11.26 Update: 文末补充C语言实现(C11标准)
向前欧拉法
{ y i + 1 = y i + h i f ( x i , y i ) y 0 = y ( a ) \left \{ \begin{array}{lr} y_{i+1}=y_i+h_i f(x_i,y_i) \\ y_0=y(a) \end{array} \right . {yi+1=yi+hif(xi,yi)y0=y(a)
from pylab import *
import warnings
warnings.filterwarnings('ignore')
求解初值问题
{
y
′
=
x
−
y
+
1
0
≤
x
≤
1
y
(
0
)
=
1
\left \{ \begin{array}{lr} y'=x-y+1 & 0\leq x \leq 1 \\ y(0)=1 \end{array} \right .
{y′=x−y+1y(0)=10≤x≤1
该微分方程的精确解为: y = x + e x p ( − x ) y=x+exp(-x) 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)
三阶龙格库塔法
{ y i + 1 = y i + 1 6 ( k 1 + 2 k 2 + k 3 ) k 1 = h f ( x i , y i ) k 2 = h f ( x i + 1 2 h , y i + 1 2 k 1 ) k 3 = h f ( x i + h , y i − k 1 + 2 k 2 ) y 0 = y ( a ) , i = 0 , 1 , ⋯ , n − 1 \left \{ \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} \right . ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧yi+1=yi+61(k1+2k2+k3)k1=hf(xi,yi)k2=hf(xi+21h,yi+21k1)k3=hf(xi+h,yi−k1+2k2)y0=y(a),i=0,1,⋯,n−1
四阶龙格库塔法
{ y i + 1 = y i + 1 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) k 1 = h f ( x i , y i ) k 2 = h f ( x i + 1 2 h , y i + 1 2 k 1 ) k 3 = h f ( x i + 1 2 h , y i + 1 2 k 2 ) k 4 = h f ( x i + h , y i + k 3 ) y 0 = y ( a ) , i = 0 , 1 , ⋯ , n − 1 \left \{ \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} \right . ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧yi+1=yi+61(k1+2k2+2k3+k4)k1=hf(xi,yi)k2=hf(xi+21h,yi+21k1)k3=hf(xi+21h,yi+21k2)k4=hf(xi+h,yi+k3)y0=y(a),i=0,1,⋯,n−1
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
改进欧拉法
{ y p = y i + h f ( x i , y i ) y i + 1 = y i + h 2 [ f ( x i , y i ) + f ( x i , y p ) ] i = 0 , 1 , ⋯ , n − 1 y 0 = y ( a ) \left \{ \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} \right . ⎩⎨⎧yp=yi+hf(xi,yi)yi+1=yi+2h[f(xi,yi)+f(xi,yp)]y0=y(a)i=0,1,⋯,n−1
求解初值问题
{
y
′
=
−
y
(
0
≤
x
≤
1
)
y
(
0
)
=
1
\left \{ \begin{array}{lr} y'=-y &(0 \leq x \leq 1) \\ y(0)=1 \\ \end{array} \right .
{y′=−yy(0)=1(0≤x≤1)
精确解为
y
=
e
x
p
(
−
x
)
y=exp(-x)
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;
}