隐式龙格库塔法举例说明

题目

用隐式中点公式求解常微分方程:
$$\{\begin{array}{l}\frac{dy}{dx}=y,\\ y(0)=1.\end{array}$$

具体分析

前期准备

首先对和在区间上进行离散化，然后构造递推公式，再进一步得到在这些位置的近似取值。
这些离散位置可以表示为:
$$x_1,x_2,\cdots ,x_n,\cdots$$
显然，$y(x)$在这些离散的点处精确取值为：
$$y(x_1),y(x_2),\cdots ,y(x_n),\cdots$$
令相邻两个离散点的间隔为一定值$h$，则有$x_n=x_1+nh(n=1,2,\cdots )$，记$f(x,y)=y$。

p级隐式R-K方法(类似于龙格-库塔法的思想)
$$\{\begin{array}{l}{y}_{n+1}=y_n+h\sum _{r=1}^p{\lambda }_r{K}_r,\\ K_r=f(x_n+{\alpha }_{r}h,y_n+h\sum _{i=1}^p{\lambda }_{ri}{K}_i),r=1,2,\cdots ,p.\end{array}$$
我们可利用$Taylor$展开确定系数。(下面以r=1为例)

确定系数

$$\{\begin{array}{ll}{y}_{n+1}=y_n+h\lambda K,& (1)\\ K=f(x_n+\alpha h,y_n+h{\lambda }_{\ast }K).& (2)\end{array}$$
首先，我们对二式中对$f$在$x_n$出$Taylor$展开，得：
$$K=f(x_n,y_n)+f_x(x_n,y_n)\alpha h+f_y(x_n,y_n)h{\lambda }_{\ast }K+o(h^2)$$
解出$K=\frac{f(x_n,y_n)+f_x(x_n,y_n)\alpha h}{1-f_y(x_n,y_n)h{\lambda }_{\ast }}$，
代入(1)式得：
$$y_{n+1}=y_n+h\lambda \frac{f(x_n,y_n)+f_x(x_n,y_n)\alpha h}{1-f_y(x_n,y_n)h{\lambda }_{\ast }}$$
与$y(x_{n+1})$在$x_n$处$Taylor$展开进行比较:
$$y(x_{n+1})=y(x_n)+\frac{h}{1!}{y}^{{}^{\prime }}(x_n)+\frac{h^2}{2!}{y}^{{}^{\prime \prime }}(x_n)+\cdots$$
得到：
$$h\lambda \frac{f(x_n,y_n)+f_x(x_n,y_n)\alpha h}{1-f_y(x_n,y_n)h{\lambda }_{\ast }}\approx \frac{h}{1!}{y}^{{}^{\prime }}(x_n)+\frac{h^2}{2!}{y}^{{}^{\prime \prime }}(x_n)$$
即：
$$\lambda \frac{f(x_n,y_n)+f_x(x_n,y_n)\alpha h}{1-f_y(x_n,y_n)h{\lambda }_{\ast }}\approx f(x_n,y_n)+\frac{h}{2}[f_x(x_n,y_n)+f_y(x_n,y_n)f(x_n,y_n)]$$
简记为：
$$\lambda \frac{f+f_x\alpha h}{1-f_{y}h{\lambda }_{\ast }}\approx f+\frac{h}{2}(f_x+f_{y}f)$$
化简得：
$$\lambda f+\lambda f_x\alpha h=f-f{f}_{y}h{\lambda }_{\ast }+\frac{1}{2}h{f}_x+\frac{1}{2}h{f}_{y}f$$
解得：
$$\lambda =1,{\lambda }_{\ast }=\frac{1}{2},\alpha =\frac{1}{2}.$$
从而得到隐式中点公式：
$$\{\begin{array}{ll}{y}_{n+1}=y_n+hK,& (1)\\ K=f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}hK).& (2)\end{array}$$

MATLAB求解

% Implicit kutta
clear all,close all,clc
f=@(x,y)y;
N=10;h=1/N;
x=linspace(0,1,N+1);
y=[1,zeros(1,N)];
K=randn(1);
for n=1:N
    K=f(x(n)+1/2*h,y(n)+1/2*h*K);
    y(n+1)=y(n)+h*K;
end
y_exact=@(x)exp(x);
plot(x,y_exact(x),'-b',x,y,'*r')
axis([0 1 0.9 3.5]),xlabel x,ylabel y
legend('Exact','ImKutta')
error=y_exact(x)-y;