隐式 Runge-Kutta (Gauss 2s)

本文用 Gauss 2 级 4 阶方法（属于隐式 Runge-Kutta 法）求解ODE，需要用到方程组的牛顿迭代，然后用 MATLAB 编程，最后输出结果的误差阶为 4。

Runge-Kutta 方法

对于常微分方程

$$y^{\prime}=f(x, y), \quad y(a)=\eta, \quad f : \mathbb{R} \times \mathbb{R}^{m} \rightarrow \mathbb{R}^{m}.$$

一般 s 级 Runge-Kutta 方法：

$$(\rm I)\left\{\begin{aligned} &y_{n+1}=y_{n}+h \sum_{i=1}^{s} b_{i} k_{i}, \\ &k_{i}=f(x_{n}+c_{i} h, y_{n}+h \sum_{j=1}^{s} a_{i j} k_{j}), \quad i=1,2, \ldots, s. \end{aligned}\right.$$

对应的 Butcher 阵：

$$\begin{array}{c|cccc} {c_{1}} & {a_{11}} & {a_{12}} & {\cdots} & {a_{1 s}} \\ {c_{2}} & {a_{21}} & {a_{22}} & {\cdots} & {a_{2 s}} \\ {\vdots}& {\vdots} & {} & {} & {\vdots} \\ {c_{s}} & {a_{s 1}} & {a_{s 2}} & {\cdots} & {a_{s s}} \\ \hline {} & b_{1} & b_{2} & \cdots & b_{s}\\ \end{array}$$ $$c=\left[c_{1}, c_{2}, \ldots, c_{s}\right]^{\mathrm{T}}, \quad b=\left[b_{1}, b_{2}, \ldots, b_{s}\right]^{\mathrm{T}}, \quad A=\left(a_{i j}\right)_{s}.$$ $$(\rm II)\left\{\begin{aligned} &y_{n+1}=y_{n}+h \sum_{i=1}^{s} b_{i} f\left(x_{n}+c_{i} h, Y_{i}\right), \\ &Y_{i}=y_{n}+h \sum_{j=1}^{s} a_{i j} f\left(x_{n}+c_{j} h, Y_{j}\right), \quad i=1,2, \ldots, s. \end{aligned}\right.$$

形式 (I) 与形式 (II) 等价，可以通过变换 得到。

Gauss Method 2 级 4 阶

$$\left\{ \begin{aligned} &y_{n+1}=y_{n}+h \sum_{i=1}^{2} b_{i} f\left(x_{n}+c_{i} h, Y_{i}\right), \\ &Y_{1}=y_{n}+h \sum_{j=1}^{2} a_{1 j} f\left(x_{n}+c_{j} h, Y_{j}\right), \\ &Y_{2}=y_{n}+h \sum_{j=1}^{2} a_{2 j} f\left(x_{n}+c_{j} h, Y_{j}\right), \end{aligned}\right.$$

对应的 Butcher 阵：

$$\begin{array}{c|cc} \frac{1}{2}-\frac{\sqrt{3}}{6} & \frac{1}{4} & \frac{1}{4}-\frac{\sqrt{3}}{6} \\[6pt] \frac{1}{2}+\frac{\sqrt{3}}{6} & \frac{1}{4}+\frac{\sqrt{3}}{6} & \frac{1}{4} \\[6pt] \hline & \frac{1}{2} & \frac{1}{2} \end{array}$$

方程组的牛顿迭代

$$\left\{\begin{aligned} &f_{1}\left(x_{1}, x_{2}, \cdots, x_{n}\right)=0, \\ &f_{2}\left(x_{1}, x_{2}, \cdots, x_{n}\right)=0, \\ &\quad \vdots \\ &f_{n}\left(x_{1}, x_{2}, \cdots, x_{n}\right)=0, \end{aligned}\right.$$

记 ，，方程组可写为

$$\boldsymbol{F}(\boldsymbol{x})=\mathbf{0}.$$

方程组的牛顿迭代法

$$\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}-\boldsymbol{F}^{\prime}(\boldsymbol{x}^{(k)})^{-1} \boldsymbol{F}\left(\boldsymbol{x}^{(k)}\right), \quad k=0,1, \cdots，$$

这里 是 Jacobi 矩阵。

ODE 例子

$$\left\{\begin{aligned} &\frac{\mathrm{d}u}{\mathrm{d}t}=f(t,u),\quad 0 < t \leqslant 1, \\ &u(0)=1. \end{aligned} \right.$$

其中：.

方程的真解：.

计算误差阶

% IRK2s_error.m
% Implicit Runge-Kutta(Gauss method) 2 stage and order 4
% u'=u in [0,1] with initial condition u(0)=1
% exact solution: ue=exp(x)
clear all; close all;
Nvec=[10 50 100 200 500 1000];
Error=[];
for n=1:length(Nvec)
    N=Nvec(n);
    h=1/N;
    x=[0:h:1];
    u(1)=1;
    Y=[1;1];
    % Newton iteration
    for i=1:N
        k=u(i);  tol=1;
        while tol>1.0e-10
            X=Y;
            D=[1-0.25*h,-h*(0.25-(sqrt(3))/6);...
            -h*( 0.25+(sqrt(3))/6),1-h*0.25];    % Jacobian matrix
            F=[X(1)-k-h*(0.25*X(1)+(0.25-(sqrt(3))/6)*X(2));...
            X(2)-k-h*((0.25+(sqrt(3))/6)*X(1)+0.25*X(2))];   % RHS
            Y=X-D\F;
            tol=norm(Y-X);
        end
        u(i+1)=k+(h/2)*(Y(1)+Y(2));
    end
    ue=exp(x);              % exact solution
    error=max(abs(u-ue));   % maximum error
    Error=[Error,error];
end
plot(log10(Nvec),log10(Error),'ro-','MarkerFaceColor','w','LineWidth',1)
hold on,
plot(log10(Nvec), log10(Nvec.^(-4)), '--')
grid on,
% title('Convergence order of Gauss method ','fontsize',12)
set(gca,'fontsize',12)
xlabel('log_{10}N','fontsize',14), ylabel('log_{10}Error','fontsize',14)

% add annotation of slope
ax = [0.62 0.58];
ay = [0.72 0.66];
annotation('textarrow',ax,ay,'String','slope = -4 ','fontsize',14)

% computating convergence order
for i=1:length(Nvec)-1     % computating convergence order
    order(i)=-log(Error(i)/Error(i+1))/(log(Nvec(i)/Nvec(i+1)));
end
Error
order

% print -dpng -r600 IRK2s_error.png
% print -depsc2 IRK2s_error.eps

输出结果

Error =
   1.0e-06 *
    0.3778    0.0006    0.0000    0.0000    0.0000    0.0000

order =
    4.0004    4.0001    3.9990    3.9310    4.3726

另一个例子

$$\left\{\begin{aligned} &\frac{\mathrm{d}u}{\mathrm{d}t}=f(t,u),\quad 0 < t \leqslant 1, \\ &u(0)=0. \end{aligned} \right.$$

其中：.

方程的真解：.

% IRK2s2_error.m
% Implicit Runge-Kutta(Gauss method) 2 stage and order 4
% u'=u^2/2-1/2 in [0,1] with initial condition u(0)=0
% exact solution: ue=(1-exp(x))/(1+exp(x))
clear all; close all;
Nvec=[10 50 100 200 500];
Error=[];
for n=1:length(Nvec)
    N=Nvec(n);
    h=1/N;
    x=[0:h:1];
    u(1)=0;
    Y=[1;1];
    % Newton iteration
    for i=1:N
        k=u(i);  tol=1;
        while tol>1.0e-10
            X=Y;
            D=[1-h*0.25*X(1), -h*(0.25-(sqrt(3))/6)*X(2);...
                -h*(0.25+(sqrt(3))/6)*X(1), 1-h*0.25*X(2)];
            F=[X(1)-k-h*(0.25*(0.5*(X(1))^2-0.5)+(0.25-(sqrt(3))/6)*(0.5*(X(2))^2-0.5));...
                X(2)-k-h*((0.25+(sqrt(3))/6)*(0.5*(X(1))^2-0.5)+0.25*(0.5*(X(2))^2-0.5))];
            Y=X-D\F;
            tol=norm(Y-X);
        end
        u(i+1)=k+(h/2)*(0.5*Y(1)^2-0.5+0.5*Y(2)^2-0.5);
    end
    ue=(1-exp(x))./(1+exp(x));  % exact solution
    error=max(abs(u-ue));       % maximum error
    Error=[Error,error];
end
plot(log10(Nvec),log10(Error),'ro-','MarkerFaceColor','w','LineWidth',1)
hold on,
plot(log10(Nvec),log10(Nvec.^(-4)),'--')
grid on,
% title('Convergence order of Gauss method ','fontsize',12)
set(gca,'fontsize',12)
xlabel('log_{10}N','fontsize',14), ylabel('log_{10}Error','fontsize',14)

% add annotation of slope
ax = [0.62 0.58];
ay = [0.72 0.66];
annotation('textarrow',ax,ay,'String','slope = -4 ','fontsize',14)

% computating convergence order
for i=1:length(Nvec)-1     % computating convergence order
    order(i)=-log(Error(i)/Error(i+1))/(log(Nvec(i)/Nvec(i+1)));
end
Error
order

% print -dpng -r600 IRK2s2_error.png
% print -depsc2 IRK2s2_error.eps

输出结果

Error =
   1.0e-08 *
    0.8913    0.0014    0.0001    0.0000    0.0000

order =
    3.9999    3.9999    4.0011    3.6619