[數學]龍格-庫塔法

原理思想

要想求出非常近似的值，有種神器叫做泰勒公式 。泰勒給出了任意一個函數都可以用多項式逼近的方法求出函數值。這與常微分方程的數值方法的思想類似，就是已知初始值，借助導數這個工具，將其近似成求另一個點的坐標。龍格-庫塔方法是用斜率的權重最后得到一個較好的斜率，然后求出函數值。

公式推導

以二階推導為例
由泰勒公式：$$f(x)=\frac{f(x_0)}{0!}+\frac{f^{'}(x_0)}{1!}(x-x_0)+\frac{f^{''}(x_0)}{2!}(x-x_0)^2+...+\frac{f^{n}(x_0)}{n!}(x-x0)^n+R_n(x) \tag{1}$$
令\(x-x_0=h\),於是得\(x=h+x_0\),(1)式即可寫成$$f(x_0+h)=\frac{f(x_0)}{0!}+hf^{'}(x_0)+h^2\frac{f^{''}(x_0)}{2!}+...+h^n\frac{f^{n}(x_0)}{n!}+R_n(h+x_0) \tag{2}$$
根據實際含義，已知初始點\((x_0,y_0)\),導函數\(y^{'}=f(x,y)\),則下個點\(y_{i+1}\)為多少？另h為步長，\(h=x_{i+1}-x_i\),所以(2)式,忽略到高階項，可以寫成$$y_{i+1}=y_i+hf(x_i,y_i)+\frac{h^2}{2!}f^{'}(x_i,y_i)+\mathcal{O}(h^3) \tag{}$$
因為y是x的函數，對x再求導，得到$$f^{'}(x_i,y_i)=\frac{\partial{f(x_i,y_i)}} {\partial{x}} + \frac{\partial{f(x_i,y_i)}} {\partial{y}}\frac{\mathrm{d}y}{\mathrm{d}x}$$
現在考慮二元的泰勒公式，有$$F(x+h,y+k)=F(x,y)+h\frac{\partial{F}}{\partial{x}}+k\frac{\partial{F}}{\partial{y}}+\mathcal{O}{(h^2,k^2)} \tag{3}$$
設\(k2=f(x_i+ph,y_i+qk_1h)\),則根據(3)式，得到$$k2=f(x_i+ph,y_i+qk_1h)=f(x_i,y_i)+ph\frac{\partial{f}}{\partial{x}}+qk_1h\frac{\partial{f}}{\partial{y}}+\mathcal{O}(h^2)$$
根據龍格-庫塔的思想(斜率權重作為最終的斜率)，有\(y_{i+1}=y_i+h(ak_1+bk_2)\)，其中a，b是斜率的權重，帶入\(k1=f(x_i,y_i)\)和k2,有$$y_{i+1}=y_i+haf(x_i,y_i)+hb\big[f(x_i,y_i)+ph\frac{\partial{f}}{\partial{x}}+qf(f_i,y_i)h\frac{\partial{f}}{\partial{y}}+\mathcal{O}(h^2) \big]$$
整理得:$$y_{i+1}=y_i+(a+b)hf(x_i,y_i)+h^2\big(bp\frac{\partial{f}}{\partial{x}}+bq\frac{\partial{f}}{\partial{y}}f(x_i,y_i)\big) + \mathcal{O}(h^3) \tag{}$$
比較這兩個式子，()和(**),系數相等得到，$$a+b=1\bp=1/2\bq=1/2$$
這樣四個未知數，三個方程，不難得到無窮多解，選擇b=1/2,so a=1/2 p=q=1,於是$$y_{i+1}=y_i+h(\frac{1}{2}k_1+\frac{1}{2}k_2)$$
高階的龍格庫塔也是如此，先用泰勒公式求出三階的或者四階方程，然后求出每個斜率的表達式，最后比較系數得出方程組

常見的四階龍格-庫塔公式：$$k_1=f(x_i,y_i)\k_2=f(x_i+0.5h,y_i+0.5k_1 h)\k_3=f(x_i+0.5h,y_i+0.5k_2 h)\k_4=f(x_i+h,y_i+k_3 h)\y_{i+1}=y_i+h\big(\frac{1}{6}k_1+\frac{1}{3}k_2+\frac{1}{3}k_3+\frac{1}{6}k_4\big)$$

MATLAB 代碼

fun = @(x,y) (x+y);

myans = rk_method(fun,0,2,1,0.25);
fprintf('result is %f\n',myans);

hold on;

% 准確值
xx = 0:0.25:1;
yy = 3*exp(xx)-xx-1;
plot(xx,yy,'r');

legend('近似值','point','准確值');
function re = rk_method(f,x0,y0,x,h)
    n = round((x-x0)/h);
    x = x0;
    y = y0;
    
    xx = [];
    yy = [];
    xx(1) = x;
    yy(1) = y;
    for i=1:n
        k1 = f(x,y);
        k2 = f(x+0.5*h,y+0.5*k1*h);
        k3 = f(x+0.5*h,y+0.5*k2*h);
        k4 = f(x+h,y+k3*h);
        
        y  =  y+h*((1/6)*k1+(1/3)*k2+(1/3)*k3+(1/6)*k4);
        x = x+h;
        xx(i+1)=x;
        yy(i+1)=y;
    end
    
    re = y;
    plot(xx,yy,'b');
    hold on;
    scatter(xx,yy,'b*');

end

Result:

可以看出和matlab自帶的函數圖像已經重合了