MATLAB常微分方程数值解

作者：凯鲁嘎吉 - 博客园 http://www.cnblogs.com/kailugaji/

1.一阶常微分方程初值问题

2.欧拉法

3.改进的欧拉法

4.四阶龙格库塔方法

5.例题

用欧拉法，改进的欧拉法及4阶经典Runge-Kutta方法在不同步长下计算初值问题。步长分别为0.2,0.4,1.0.

matlab程序：

1function z=f(x,y)
2z=-y*(1+x*y);
3
4function R_K(h)
5%欧拉法
6y=1;
7fprintf('欧拉法：x=%f, y=%f\n',0,1);
8for i=1:1/h
9    x=(i-1)*h;
10    K=f(x,y);
11    y=y+h*K;
12    fprintf('欧拉法：x=%f, y=%f\n',x+h,y);
13end
14fprintf('\n');
15%改进的欧拉法
16y=1;
17fprintf('改进的欧拉法：x=%f, y=%f\n',0,1);
18for i=1:1/h
19    x=(i-1)*h;
20    K1=f(x,y);
21    K2=f(x+h,y+h*K1);
22    y=y+(h/2)*(K1+K2);
23    fprintf('改进的欧拉法：x=%f, y=%f\n',x+h,y);
24end
25fprintf('\n');
26 %龙格库塔方法
27 y=1;
28fprintf('龙格库塔法：x=%f, y=%f\n',0,1);
29for i=1:1/h
30    x=(i-1)*h;
31    K1=f(x,y);
32    K2=f(x+h/2,y+(h/2)*K1);
33    K3=f(x+h/2,y+(h/2)*K2);
34    K4=f(x+h,y+h*K3);
35    y=y+(h/6)*(K1+2*K2+2*K3+K4);
36    fprintf('龙格库塔法：x=%f, y=%f\n',x+h,y);
37end

结果：

1>> R_K(0.2)
2欧拉法：x=0.000000, y=1.000000
3欧拉法：x=0.200000, y=0.800000
4欧拉法：x=0.400000, y=0.614400
5欧拉法：x=0.600000, y=0.461321
6欧拉法：x=0.800000, y=0.343519
7欧拉法：x=1.000000, y=0.255934
8
9改进的欧拉法：x=0.000000, y=1.000000
10改进的欧拉法：x=0.200000, y=0.807200
11改进的欧拉法：x=0.400000, y=0.636118
12改进的欧拉法：x=0.600000, y=0.495044
13改进的欧拉法：x=0.800000, y=0.383419
14改进的欧拉法：x=1.000000, y=0.296974
15
16龙格库塔法：x=0.000000, y=1.000000
17龙格库塔法：x=0.200000, y=0.804636
18龙格库塔法：x=0.400000, y=0.631465
19龙格库塔法：x=0.600000, y=0.489198
20龙格库塔法：x=0.800000, y=0.377225
21龙格库塔法：x=1.000000, y=0.291009
22>> R_K(0.4)
23欧拉法：x=0.000000, y=1.000000
24欧拉法：x=0.400000, y=0.600000
25欧拉法：x=0.800000, y=0.302400
26
27改进的欧拉法：x=0.000000, y=1.000000
28改进的欧拉法：x=0.400000, y=0.651200
29改进的欧拉法：x=0.800000, y=0.405782
30
31龙格库塔法：x=0.000000, y=1.000000
32龙格库塔法：x=0.400000, y=0.631625
33龙格库塔法：x=0.800000, y=0.377556
34>> R_K(1)
35欧拉法：x=0.000000, y=1.000000
36欧拉法：x=1.000000, y=0.000000
37
38改进的欧拉法：x=0.000000, y=1.000000
39改进的欧拉法：x=1.000000, y=0.500000
40
41龙格库塔法：x=0.000000, y=1.000000
42龙格库塔法：x=1.000000, y=0.303395

注意：在步长h为0.4时，要将for i=1:1/h改为for i=1:0.8/h。