MATLAB常微分方程的數值解法
作者：凱魯嘎吉 - 博客園
http://www.cnblogs.com/kailugaji/
一、實驗目的
科學技術中常常要求解常微分方程的定解問題，所謂數值解法就是求未知函數在一系列離散點處的近似值。
二、實驗原理
三、實驗程序
1. 尤拉公式程序
四、實驗內容
選一可求解的常微分方程的定解問題，分別用以上1, 4兩種方法求出未知函數在
節點處的近似值，並對所求結果與分析解的（數值或圖形）結果進行比較。
五、解答
1. 程序
求解初值問題
取n=10
源程序：
euler23.m:
function [A1,A2,B1,B2,C1,C2]=euler23(a,b,n,y0)
%歐拉法解一階常微分方程
%初始條件y0
h = (b-a)/n; %步長h
%區域的左邊界a
%區域的右邊界b
x = a:h:b;
m=length(x);
%前向歐拉法
y = y0;
for i=2:m
y(i)=y(i-1)+h*oula(x(i-1),y(i-1));
A1(i)=x(i);
A2(i)=y(i);
end
plot(x,y,'r-');
hold on;
%改進歐拉法
y = y0;
for i=2:m
y(i)=y(i-1)+h/2*( oula(x(i-1),y(i-1))+oula(x(i),y(i-1))+h*(oula(x(i-1),x(i-1))));
B1(i)=x(i);
B2(i)=y(i);
end
plot(x,y,'m-');
hold on;
%歐拉兩步公式
y=y0;
y(2)=y(1)+h*oula(x(1),y(1));
for i=2:m-1
y(i+1)=y(i-1)+2*h*oula(x(i),y(i));
C1(i)=x(i);
C2(i)=y(i);
end
plot(x,y,'b-');
hold on;
%精確解用作圖
xx = x;
f = dsolve('Dy=-3*y+8*x-7','y(0)=1','x');%求出解析解
y = subs(f,xx); %將xx代入解析解，得到解析解對應的數值
plot(xx,y,'k--');
legend('前向歐拉法','改進歐拉法','歐拉兩步法','解析解');
oula.m:
function f=oula(x,y)
f=-3*y+8*x-7;
2. 運算結果
A1,A2為前向歐拉法在節點處的近似值，B1,B2為改進的歐拉法在節點處的近似值，C1,C2為歐拉公式法在節點處的近似值。
>> [A1,A2,B1,B2,C1,C2]=euler23(0,1,10,1)
A1 =
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000
A2 =
0 0 -0.6200 -0.9740 -1.1418 -1.1793 -1.1255 -1.0078 -0.8455 -0.6518 -0.4363
B1 =
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000
B2 =
0 0.0050 -0.6090 -0.9563 -1.1169 -1.1468 -1.0853 -0.9597 -0.7893 -0.5875 -0.3638
C1 =
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000
C2 =
0 0 -0.2400 -0.9360 -0.5984 -1.3370 -0.3962 -1.5392 0.2473 -1.8076
>> [A1,A2,B1,B2,C1,C2]=euler23(0,1,10,1)
A1 =
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000
A2 =
0 0 -0.6200 -0.9740 -1.1418 -1.1793 -1.1255 -1.0078 -0.8455 -0.6518 -0.4363
B1 =
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000
B2 =
0 0.0050 -0.6090 -0.9563 -1.1169 -1.1468 -1.0853 -0.9597 -0.7893 -0.5875 -0.3638
C1 =
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000
C2 =
0 0 -0.2400 -0.9360 -0.5984 -1.3370 -0.3962 -1.5392 0.2473 -1.8076
3. 拓展（方法改進、體會等）
從以上圖形可以看出，在n=10時，改進的歐拉法精度更高，而歐拉兩步法所求結果震盪不收斂，越接近1，震盪幅度越大，於是取n=100,時，結果如下所示：
當n=1000時，結果如下圖：
當n=100時，三種方法與解析解非常接近，當n=1000時，幾乎四者位於一條線中，從實驗結果看出，n越大時，結果越精確。