Andy123t
0%

隐式 Runge-Kutta (Gauss 2s)

发表于 更新于 分类于

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

Runge-Kutta 方法

对于常微分方程

一般 s 级 Runge-Kutta 方法：

对应的 Butcher 阵：

形式 (I) 与形式 (II) 等价，可以通过变换 $ k_{i}=f\left(x_{n}+c_{i} h, Y_{i}\right),~i=1,2, \dots, s$ 得到。

Gauss Method 2 级 4 阶

对应的 Butcher 阵：

方程组的牛顿迭代

记 $\boldsymbol{x}=\left(x_{1}, x_{2}, \cdots, x_{n}\right)^{\mathrm{T}} \in \mathbb{R}^{n}$，$\boldsymbol{F}=\left(f_{1}, f_{2}, \cdots, f_{n}\right)^{\mathrm{T}}$，方程组可写为

方程组的牛顿迭代法

这里 $\boldsymbol{F}^{\prime}(\boldsymbol{x})$ 是 Jacobi 矩阵。

ODE 例子

其中：$f(t,u)=u$.

方程的真解：$u(x)=e^{x}$.

计算误差阶

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
% 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

输出结果

1
2
3
4
5
6
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

IRK2s_error.png

另一个例子

其中：$f(t,u)=\displaystyle \frac{u^2}{2}+\frac{1}{2}$.

方程的真解：$u(x)=\displaystyle \frac{1-e^{x}}{1+e^{x}}$.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
% 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

输出结果

1
2
3
4
5
6
Error =
   1.0e-08 *
    0.8913    0.0014    0.0001    0.0000    0.0000

order =
    3.9999    3.9999    4.0011    3.6619

IRK2s2_error.png