matlab解常微分方程常用数值解法2：龙格库塔方法_用标准四阶龙格-库塔方法求解常微分方程初值问题matlab程序实现

总结和记录一下matlab求解常微分方程常用的数值解法，本文将介绍龙格库塔方法（Runge-Kutta Method）。

龙格库塔迭代的基本思想是：

$$x_{k+1}=x_k+a{k}_1+b{k}_2$$

$$k_1=hf\left(x_k,t_k\right)\text{ and }{k}_2=hf\left(x_k+B_1{k}_1,t_k+A_{1}h\right)$$

其中$a,b,A_1,B_1$是未知的，我们进行推导。

首先对$f(x+k,t+h)$进行泰勒展开：

f ( x + k , t + h ) = f ( x , t ) + ( k f x + h f t ) + 1 2 ( k 2 f x x + 2 k h f x t + h 2 f t t ) + 1 6 ( k 3 f x x x + 3 k 2 h f x x t + 3 k h 2 f x t t + h 3 f t t t ) + ⋯

$$\begin{aligned} f(x+k, t+h) & =f(x, t)+\left(k f_{x}+h f_{t}\right)+\frac{1}{2}\left(k^{2} f_{xx}+2 k h f_{xt }+h^{2} f_{tt}\right) \\ & +\frac{1}{6}\left(k^{3} f_{xxx}+3 k^{2} h f_{xxt}+3 k h^{2} f_{xtt}+h^{3} f_{ttt}\right)+\cdots \end{aligned}$$ f(x+k,t+h)​=f(x,t)+(kfx​+hft​)+21​(k2fxx​+2khfxt​+h2ftt​)+61​(k3fxxx​+3k2hfxxt​+3kh2fxtt​+h3fttt​)+⋯​
利用泰勒展开我们展开$k_2$：
$$\begin{array}{rl}{k}_2& =hf\left[x_k+B_{1}hf\left(x_k,t_k\right),t_k+A_{1}h\right]\\ & =h\left[f\left(x_k,t_k\right)+\left(B_{1}hf{f}_x+A_{1}h{f}_t\right)\right]\\ & =hf+A_1{h}^2{f}_t+B_1{h}^{2}f{f}_x,\end{array}$$
上式中的$f$为$f(x_k,t_k)$，我们将上式代回到最开始的式子$x_{k+1}=x_k+a{k}_1+b{k}_2$得到(*)式

$$x_{k+1}=x_k+(a+b)hf+\left(A_{1}b{f}_t+B_{1}bf{f}_x\right)h^2(\ast )$$

对应于二阶泰勒展式：

$$x_{k+1}=x_k+h{x}_k^{\prime }+\frac{1}{2}{h}^2{x}_k^{\prime \prime }$$

对于微分方程我们知道$x^{\prime }=f(x,t)$，于是对$t$求二阶导数：

$$x^{\prime \prime }=f_t+f_x{x}^{\prime }=f_t+f{f}_x$$

于是有：

$$x_{k+1}=x_k+hf+\frac{1}{2}{h}^2\left(f_t+f{f}_x\right)$$

对比(*)式我们有：

$$a+b=1,A_{1}b=\frac{1}{2},\text{ and }{B}_{1}b=\frac{1}{2}\text{. }$$

如果我们取$a=b=\frac{1}{2}$，那么$A_1=B_1=1$，为二阶的龙哥库塔算法。其形式和改进的欧拉法一致:

$$x_{k+1}=x_k+\frac{1}{2}\left(k_1+k_2\right)$$

从二阶的出发，我们可以得到改进的一套更好的框架，一个常用的龙哥库塔方法是四阶的：

$$x_{k+1}=x_k+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)\text{, }$$

其中：

k 1 = h f ( x k , t k ) k 2 = h f ( x k + 1 2 k 1 , t k + 1 2 h ) k 3 = h f ( x k + 1 2 k 2 , t k + 1 2 h ) k 4 = h f ( x k + k 3 , t k + h )

$$\begin{aligned} &k_{1}=h f\left(x_{k},t_{k} \right) \\ &k_{2}=h f\left(x_{k}+\frac{1}{2} k_{1},t_{k}+\frac{1}{2} h \right) \\ &k_{3}=h f\left(x_{k}+\frac{1}{2} k_{2},t_{k}+\frac{1}{2} h \right) \\ &k_{4}=h f\left(x_{k}+k_{3},t_{k}+h\right) \end{aligned}$$ ​k1​=hf(xk​,tk​)k2​=hf(xk​+21​k1​,tk​+21​h)k3​=hf(xk​+21​k2​,tk​+21​h)k4​=hf(xk​+k3​,tk​+h)​

例子

$$x^{\prime }=x+t,x(0)=1$$

clc
clear
% 测试4个不同的步长
test_times = 3;
h_test = [0.10, 0.05, 0.01];%根据步长数改变

% 保存时间、差分时间和步长
h_res=ones(1,test_times);
t_res=cell(1,test_times);%时间
tplot_res=cell(1,test_times);%差分的时间，比时间长度少1
% 保存两种数值方法和解析解的计算结果
x_rk_res=cell(1,test_times);
x_exact_res=cell(1,test_times);
% 保存误差
diff_res=cell(1,test_times);



for i = 1:test_times
% 设置步长间隔和步长数
    h = h_test(i);
    n = 10/h;
% 设置初始条件
    t=zeros(n+1,1); t(1) = 0;
    x_rk=zeros(n+1,1); x_rk(1) = 1;
    x_exact=zeros(n+1,1); x_exact(1) = 1;
% 设置龙哥库塔方法误差存放向量（和精确解比较）
    diff = zeros(n,1); tplot = zeros(n,1);
% 定义微分方程
    f = @(tt,xx)(xx+tt);
    for k = 1:n
        x_local = x_rk(k); t_local = t(k);
        k1 = h * f(t_local,x_local);
        k2 = h * f(t_local + h/2,x_local + k1/2);
        k3 = h * f(t_local + h/2,x_local + k2/2);
        k4 = h * f(t_local + h,x_local + k3);
        t(k+1) = t_local + h;
        x_rk(k+1) = x_local + (k1+2*k2+2*k3+k4) / 6;
        x_exact(k+1) = 2*exp(t(k+1)) - t(k+1) - 1;
        tplot(k) = t(k);
        diff(k) = x_rk(k+1) - x_exact(k+1);
        diff(k) = abs(diff(k) / x_exact(k+1));
    end
    diff_res{i}=diff;
	tplot_res{i}=tplot;
	h_res(i)=h;
	x_rk_res{i}=x_rk;
	x_exact_res{i}=x_exact;
	t_res{i}=t; 
end

figure
% 不同步长下的解析解和龙哥库塔法的结果
for i=1:test_times
    subplot(2,2,i)
    plot(t_res{i},x_exact_res{i},'r-',t_res{i},x_rk_res{i},'b--')
    xlabel('t','Fontsize',18)
    ylabel('x','Fontsize',18)
    legend({'Analytical method','Runge-Kutta method'},'Location','best')
    legend boxoff;
    title(['h = ',num2str(h_res(i))]);
    axis tight
end

% 相对误差图
subplot(2,2,4)
for i = 1:test_times
    semilogy(tplot_res{i},diff_res{i},'b-')
    hold on
    num1 = 2; num2 = 2/h_test(i);
    text(num1,diff_res{i}(num2),['h = ',num2str(h_res(i))],'Fontsize',15,...
    'HorizontalAlignment','right',...
    'VerticalAlignment','bottom')
end
xlabel('t','Fontsize',20)
ylabel('|relative error|','Fontsize',20)
title('Runge-Kutta method''s relative error')
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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

可以看到使用龙格库塔法，步长$h=0.1$精度就已经很高了。