老规矩,数学原理什么的就不写了。
直接贴代码和实例演示,以下代码基于python和numpy。
在这里,我将用代码实现改进的欧拉公式、标准(经典)的四阶龙格-库塔法(改进的 Euler 公式和标准(经典)的四阶 Runge-Kutta 法)和线性常微分方程第一边值问题的差分解法。
没时间看的同学,可以直接跳到总结。
文章目录
- 改进的 Euler 方法
- 定义函数
- 参数说明
- 实例运行
- 标准(经典)的四阶 Runge-Kutta 法
- 定义函数
- 参数说明
- 实例运行
- 线性常微分方程第一边值问题的差分解法
- 定义函数
- 参数说明
- 实例运行
- 总结
改进的 Euler 方法
什么是改进的 Euler 公式,它长什么样,如下图所示:
首先
import numpy as np
定义函数
以下便是我定义的函数:
def euler_pro(x0,f,y0,h):
x=np.arange(x0[0],x0[1]+h,h)
y=np.empty(x.shape)
y[0]=y0
for i in range(int((x0[1]-x0[0])/h)):
y[i+1]=y[i]+h/2*(f(x[i],y[i])+f(x[i+1],y[i]+h*f(x[i],y[i])))
z=np.hstack([x.reshape(-1,1),y.reshape(-1,1)])
print(z[1:])
return y[1:]
参数说明
“x0”指的是自变量的定义域。
“f”可以是y的导数与x、y的函数关系,此时y的导数在等号的左侧,x、y在等号的右侧。
“f”代表的函数就为等号的右侧。讲的不是太清楚,大致理解一下,待会儿看实例运行时就懂了。
“y0”代表初值条件,待会儿看实例运行时就知道是什么了。
“h”指的是步长。
实例运行
下图为给出的实例,我们一一解释其中我们要的参数。
函数中的“x0”就是图中的“0<x<=1”。
函数中的“f”就是图中的“y-2×x/y”。
函数中的“y0”就是图中的“y(0)=1”。
函数中的“h”就是图中的“步长h=0.1”。
实际写入操作如下。
f1=lambda x,y: y-2*x/y
x=np.array([0,1])
y0=1
h=0.1
euler_pro(x,f1,y0,h)
得出以下结果:
[[0.1 1.09590909]
[0.2 1.18409657]
[0.3 1.26620136]
[0.4 1.34336015]
[0.5 1.41640193]
[0.6 1.4859556 ]
[0.7 1.55251409]
[0.8 1.61647478]
[0.9 1.67816636]
[1. 1.7378674 ]]
这个数组代表的含义是:
x=0.1时,y=1.09590909;
x=0.2时,y=1.18409657;
x=0.3时,y=1.26620136;
……
x=1时,y=1.7378674。
注意此初值问题的解析解为y=(1+2x)0.5, 从下图(因为精确度不同,所以数值有点略微不同)可以看出, 数值解 yn与解析解 y(xn) 基本上只保证小数点后两位数的准确性(h2 = 0.01)。
标准(经典)的四阶 Runge-Kutta 法
什么是标准(经典)的四阶 Runge-Kutta 法,它长什么样,如下图所示:
定义函数
以下便是我定义的函数:
def runge_kutta_4(x0,f,y0,h):
x=np.arange(x0[0],x0[1]+h,h)
y=np.empty(x.shape)
y[0]=y0
k1=k2=k3=k4=0
for i in range(int((x0[1]-x0[0])/h)):
k1=f(x[i],y[i])
k2=f(x[i]+h/2,y[i]+h*k1/2)
k3=f(x[i]+h/2,y[i]+h*k2/2)
k4=f(x[i]+h,y[i]+h*k3)
y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6
z=np.hstack([x.reshape(-1,1),y.reshape(-1,1)])
print(z[1:])
return y[1:]
参数说明
参数意义和上一个函数一样,就不再赘述。
实例运行
运行和上面一样的实例,不过步长改为0.2。
f1=lambda x,y: y-2*x/y
x=np.array([0,1])
y0=1
h=0.2
runge_kutta_4(x,f1,y0,h)
得出以下结果:
[[0.2 1.18322929]
[0.4 1.34166693]
[0.6 1.48328146]
[0.8 1.61251404]
[1. 1.73214188]]
这个数组代表的含义是:
x=0.2时,y=1.18322929;
……
x=1时,y=1.73214188。
同上理,此初值问题的解析解为y=(1+2x)0.5, 从下图可以看出, 数值解 yn与解析解 y(xn) 能保证小数点后四位数的准确性(h2 = 0.0001)。
显然在计算量大致相同的情况下, 标准的四阶 Runge-Kutta 法比改进的 Euler 方法精确度更高。
线性常微分方程第一边值问题的差分解法
什么是线性常微分方程,它长什么样,大概如下图所示:
然后我们要求解的是在定义域[a,b]中,步长为h时,y在各个点上的取值。
定义函数
以下便是我定义的函数:
def first_yjxxcwf(q,f,x0,y0,h):
a=x0[0]
b=x0[1]
c=y0[0]
d=y0[1]
n=int((b-a)/h)
x=np.arange(a,b+h,h)
an=-(2+q(x[1:-1])*h**2)
if type(an) is np.ndarray:
pass
else:
an=np.full(n-1,an)
d1=h**2*f(x[1])-c
dn=h**2*f(x[-2])-d
dh=h**2*f(x[2:-2])
if type(dh) is np.ndarray:
pass
else:
dh=np.full(n-3,dh)
dh=np.insert(dh,0,d1)
dh=np.append(dh,dn)
A=np.zeros((n-1,n-1))
for i in range(n-1):
A[i,i]=an[i]
for i in range(n-2):
A[i,i+1]=A[i+1,i]=1
z=np.hstack([x[1:-1].reshape(-1,1),np.linalg.inv(A).dot(dh).reshape(-1,1)])
print(z)
return np.linalg.inv(A).dot(dh)
参数说明
“q”指的是图中的“q(x)”。
“f”指的是图中的“f(x)”。
“x0”指的是自变量的定义域。
“y0”指的是两个初值条件。
“h”指的是步长。
实例运行
下图为给出的实例,步长h=0.1。
q=lambda x: 1
f=lambda x: x
x=[0,1]
y=[0,1]
h=0.1
first_yjxxcwf(q,f,x,y,h)
得出以下结果:
[[0.1 0.07048938]
[0.2 0.14268365]
[0.3 0.21830476]
[0.4 0.29910891]
[0.5 0.38690416]
[0.6 0.48356844]
[0.7 0.59106841]
[0.8 0.71147907]
[0.9 0.84700451]]
这个数组代表的含义和上面两个例子一样。
总结
在计算量大致相同的情况下, 标准的 Runge-Kutta 法比改进的 Euler 方法精确度更高。
三者的代码总结如下:
# 改进的 Euler 方法
def euler_pro(x0,f,y0,h):
x=np.arange(x0[0],x0[1]+h,h)
y=np.empty(x.shape)
y[0]=y0
for i in range(int((x0[1]-x0[0])/h)):
y[i+1]=y[i]+h/2*(f(x[i],y[i])+f(x[i+1],y[i]+h*f(x[i],y[i])))
z=np.hstack([x.reshape(-1,1),y.reshape(-1,1)])
print(z[1:])
return y[1:]
# 标准的四阶 Runge-Kutta 法
def runge_kutta_4(x0,f,y0,h):
x=np.arange(x0[0],x0[1]+h,h)
y=np.empty(x.shape)
y[0]=y0
k1=k2=k3=k4=0
for i in range(int((x0[1]-x0[0])/h)):
k1=f(x[i],y[i])
k2=f(x[i]+h/2,y[i]+h*k1/2)
k3=f(x[i]+h/2,y[i]+h*k2/2)
k4=f(x[i]+h,y[i]+h*k3)
y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6
z=np.hstack([x.reshape(-1,1),y.reshape(-1,1)])
print(z[1:])
return y[1:]
# 线性常微分方程第一边值问题的差分解法
def first_yjxxcwf(q,f,x0,y0,h):
a=x0[0]
b=x0[1]
c=y0[0]
d=y0[1]
n=int((b-a)/h)
x=np.arange(a,b+h,h)
an=-(2+q(x[1:-1])*h**2)
if type(an) is np.ndarray:
pass
else:
an=np.full(n-1,an)
d1=h**2*f(x[1])-c
dn=h**2*f(x[-2])-d
dh=h**2*f(x[2:-2])
if type(dh) is np.ndarray:
pass
else:
dh=np.full(n-3,dh)
dh=np.insert(dh,0,d1)
dh=np.append(dh,dn)
A=np.zeros((n-1,n-1))
for i in range(n-1):
A[i,i]=an[i]
for i in range(n-2):
A[i,i+1]=A[i+1,i]=1
z=np.hstack([x[1:-1].reshape(-1,1),np.linalg.inv(A).dot(dh).reshape(-1,1)])
print(z)
return np.linalg.inv(A).dot(dh)
如果哪天数学作业涉及微分方程组的初值问题,那我会再编一个向量形式的标准四阶龙格-库塔。因为好像挺麻烦的样子,就不顺便编了。
如果哪天数学作业涉及二阶微分方程的数值解法,那我会再编一个差分法。