MATLAB常微分方程數值解

作者：凱魯嘎吉 - 博客園 http://www.cnblogs.com/kailugaji/

1.一階常微分方程初值問題

2.歐拉法

3.改進的歐拉法

4.四階龍格庫塔方法

5.例題

用歐拉法，改進的歐拉法及4階經典Runge-Kutta方法在不同步長下計算初值問題。步長分別為0.2,0.4,1.0.

matlab程式：

function z=f(x,y)
z=-y*(1+x*y);

function R_K(h)
%歐拉法
y=1;
fprintf('歐拉法：x=%f, y=%f\n',0,1);
for i=1:1/h
    x=(i-1)*h;
    K=f(x,y);
    y=y+h*K;
    fprintf('歐拉法：x=%f, y=%f\n',x+h,y);
end
fprintf('\n');
%改進的歐拉法
y=1;
fprintf('改進的歐拉法：x=%f, y=%f\n',0,1);
for i=1:1/h
    x=(i-1)*h;
    K1=f(x,y);
    K2=f(x+h,y+h*K1);
    y=y+(h/2)*(K1+K2);
    fprintf('改進的歐拉法：x=%f, y=%f\n',x+h,y);
end
fprintf('\n');
 %龍格庫塔方法
 y=1;
fprintf('龍格庫塔法：x=%f, y=%f\n',0,1);
for i=1:1/h
    x=(i-1)*h;
    K1=f(x,y);
    K2=f(x+h/2,y+(h/2)*K1);
    K3=f(x+h/2,y+(h/2)*K2);
    K4=f(x+h,y+h*K3);
    y=y+(h/6)*(K1+2*K2+2*K3+K4);
    fprintf('龍格庫塔法：x=%f, y=%f\n',x+h,y);
end

結果：

>> R_K(0.2)
歐拉法：x=0.000000, y=1.000000
歐拉法：x=0.200000, y=0.800000
歐拉法：x=0.400000, y=0.614400
歐拉法：x=0.600000, y=0.461321
歐拉法：x=0.800000, y=0.343519
歐拉法：x=1.000000, y=0.255934

改進的歐拉法：x=0.000000, y=1.000000
改進的歐拉法：x=0.200000, y=0.807200
改進的歐拉法：x=0.400000, y=0.636118
改進的歐拉法：x=0.600000, y=0.495044
改進的歐拉法：x=0.800000, y=0.383419
改進的歐拉法：x=1.000000, y=0.296974

龍格庫塔法：x=0.000000, y=1.000000
龍格庫塔法：x=0.200000, y=0.804636
龍格庫塔法：x=0.400000, y=0.631465
龍格庫塔法：x=0.600000, y=0.489198
龍格庫塔法：x=0.800000, y=0.377225
龍格庫塔法：x=1.000000, y=0.291009
>> R_K(0.4)
歐拉法：x=0.000000, y=1.000000
歐拉法：x=0.400000, y=0.600000
歐拉法：x=0.800000, y=0.302400

改進的歐拉法：x=0.000000, y=1.000000
改進的歐拉法：x=0.400000, y=0.651200
改進的歐拉法：x=0.800000, y=0.405782

龍格庫塔法：x=0.000000, y=1.000000
龍格庫塔法：x=0.400000, y=0.631625
龍格庫塔法：x=0.800000, y=0.377556
>> R_K(1)
歐拉法：x=0.000000, y=1.000000
歐拉法：x=1.000000, y=0.000000

改進的歐拉法：x=0.000000, y=1.000000
改進的歐拉法：x=1.000000, y=0.500000

龍格庫塔法：x=0.000000, y=1.000000
龍格庫塔法：x=1.000000, y=0.303395

註意：在步長h為0.4時，要將for i=1:1/h改為for i=1:0.8/h。