常微分方程数值解法

微分方程数值解

核心思想：用差分代替微分

常微分方程

eg1:$$\begin{equation}
\begin{cases}
y’=y-\frac{2x}{y}, 0<x<1 \\
y\left( 0 \right) =1 \
\end{cases}
\end{equation}$$

欧拉法

迭代方程：

$y_{n+1}=y_n+hf\left( x_n,y_n \right)$

将例题1微分方程转化为迭代方程：

$y_{n+1}=y_n+h\left( y_n-\frac{2x_n}{y_n} \right)$

function yn1=myfun(xn,yn,h)
yn1=yn+h*(yn-2*xn/yn);
end

直接一个循环就出来了：

clc;clear;
yn=1;%设置初值
xn=0;%设置起点
Xs=1;%设置终点
h=0.1;%设置步长
x(1)=xn;y(1)=yn;%用于存储变量与函数值
i=1;
while xn<Xs             %当xn达到研究区间上限结束
    yn=myfun(xn,yn,h);
    xn=xn+h;
    i=i+1;
    x(i)=xn;y(i)=yn;
end
plot(x,y,x,sqrt(1+2*x));
xlabel('自变量x')
ylabel('函数值y')
legend('欧拉法','实际值');
figure;
plot(x,y-sqrt(1+2*x))
legend('误差值')

结果如下：

改进欧拉法

改进欧拉法实则首先用欧拉法进行预测估计：

$\bar{y}_{n+1}=y_n+hf\left( x_n,y_n \right)$

再利用梯形公式估计：

$y_{n+1}=y_n+\frac{h}{2}\left[ f\left( x_n,y_n \right) +f\left( x_{n+1},\bar{y}_{n+1} \right) \right]$

转化为平均化形式：

$\begin{cases} y_p=y_n+hf\left( x_n,y_n \right)\\ y_c=y_n+hf\left( x_{n+1},y_p \right)\\ y_{n+1}=\frac{1}{2}\left( y_p+y_c \right)\\ \end{cases}$

相比于欧拉法只是迭代方程产生区别：

function yn1=myfun(xn,yn,h)
yp=yn+h*(yn-2*xn/yn);
yc=yn+h*(yn-2*(xn+h)/yp);
yn1=1/2*(yp+yc);
end

可以看出确实改善了求解精度

显式R-K方法

R-K方法是对欧拉方法的推广，只是将：

$\varDelta y=y_{n+1}-y_n=\int_{x_n}^{x_{n+1}}{f\left( x,y\left( x \right) \right) dx}\approx h\sum_{i=1}^r{c_if\left( x_n+\lambda _ih,y\left( x_n+\lambda _ih \right) \right)}$

积分的精度调高了，为了保证精度和减少计算，一般采用四阶龙格库塔法，通过待定系数可以求出四阶R-K方程模型：

$\begin{cases} y_{n+1}=y_n+\frac{h}{6}\left( K_1+2K_2+2K_3+K_4 \right)\\ K_1=f\left( x_n,y_n \right)\\ K_2=f\left( x_n+\frac{h}{2},y_n+\frac{h}{2}K_1 \right)\\ K_3=f\left( x_n+\frac{h}{2},y_n+\frac{h}{2}K_2 \right)\\ K_4=f\left( x_n+h,y_n+hK_3 \right)\\ \end{cases}$

注意：R-K法是基于泰勒展开的方法，要求解具有较好的光滑性质，如果解得光滑性差，使用四阶龙格库塔求得的结果效果反而不如改进的欧拉方法，建模时需要对不同问题进行判断使用

结果：

其截断误差是:

$o(h^5)$

刚性方程组

eg2：

$\begin{cases} u'=-1000.25u+999.75v+0.5\\ v'=999.75u-1000.25v+0.5\\ u\left( 0 \right) =1\\ v\left( 0 \right) =-1\\ \end{cases}$

利用龙格库塔求出结果：

我针对特殊情况总结了一下流程：

首先对方程组进行降阶处理转化为一阶n元微分方程组：

$\begin{cases} y_{1}^{'}=f_1\left( x,y_1,y_2\cdots y_N \right)\\ y_{2}^{'}=f_2\left( x,y_1,y_2\cdots y_N \right)\\ \vdots\\ y_{N}^{'}=f_N\left( x,y_1,y_2\cdots y_N \right)\\ y_1\left( x_0 \right) =a_1\\ y_2\left( x_0 \right) =a_2\\ \vdots\\ y_N\left( x_0 \right) =a_N\\ \end{cases}$

其次根据R-K公式得到：

$\begin{cases} y_{n+1}^{1}=y_{n}^{1}+\frac{h}{6}\left( K_{1}^{1}+2K_{2}^{1}+2K_{3}^{1}+K_{4}^{1} \right)\\ y_{n+1}^{2}=y_{n}^{2}+\frac{h}{6}\left( K_{1}^{2}+2K_{2}^{2}+2K_{3}^{2}+K_{4}^{2} \right)\\ \vdots\\ y_{n+1}^{N}=y_{n}^{N}+\frac{h}{6}\left( K_{1}^{N}+2K_{2}^{N}+2K_{3}^{N}+K_{4}^{N} \right)\\ \end{cases}$

其中：

$\begin{cases} K_{1}^{i}=f_i\left( x,y_1,y_2\cdots y_N \right)\\ K_{2}^{i}=f_i\left( x+\frac{h}{2},y_1+\frac{h}{2}K_{1}^{1},y_2+\frac{h}{2}K_{1}^{2},\cdots ,y_N+\frac{h}{2}K_{1}^{N} \right)\\ K_{3}^{i}=f_i\left( x+\frac{h}{2},y_1+\frac{h}{2}K_{2}^{1},y_2+\frac{h}{2}K_{2}^{2},\cdots ,y_N+\frac{h}{2}K_{2}^{N} \right)\\ K_{4}^{i}=f_i\left( x+h,y_1+hK_{3}^{1},y_2+hK_{3}^{2},\cdots ,y_N+hK_{3}^{N} \right)\\ \end{cases}$

即可