第六章 解常微分方程組(Systems of Ordinary Diff. Equations)
本章先解一階聯立常微分方程式,其數值方法與解一階常微分方程式類似. 然後再利用解一階聯立常微分方程式的方法, 解高階的常微分方程式(組). 在本章中包含Matlab 的m-file
1. taylor4sys.m 2. rk4sys.m
3. rk4sysnew.m ; fxsys.m 4. adamsmsys.m
將須要的m-file之檔案夾加入搜尋路徑中
path('d:\\numerical', path)
註: 如果你有安裝Matlab Notebook 要執行下列input cells (綠色指令敘述) 之前必須先執行上面的cell – [path (…) ] 藍色的內容是Matlab [output cells]
1. taylor4sys.m – 利用泰勒級數方法取到第4 階
X(t+h) =X(t) +h X'(t) + h2/2 X''(t) +h3/3! X'''(t) + h4/4! X(4)(t) 來估算聯立常微分方程式 X'(t)=F(X,t)的數值解
顯示taylor4sys.m的內容
type taylor4sys.m
function rs= taylor4sys(X0,a,b,m)
%to return the approximation values of X(t)
%by Taylor series method of order 4, X0 is the initial %t is in[a,b] n=length(X0); rs=zeros(n,m) ;
d1=zeros(n,1) ; d2=d1 ;d3=d1; d4=d1; t=a ; X=X0 ; h=(b-a)/m ; for k=1:m
%for your F(t,X) compute the derivatives of X', X'', X''' and X4 at t d1(1)=X(1)-X(2)+t*(2-t*(1+t)) ; % d1 is X' d1(2)=X(1)+X(2)+t^2*(-4+t) ;
d2(1)=d1(1)-d1(2)+2-t*(2+3*t) ; % d2 is X'' d2(2)=d1(1)+d1(2)+t*(-8+3*t) ;
d3(1)=d2(1)-d2(2)-2-6*t ; % d3 is X''' d3(2)=d2(1)+d2(2)-8+6*t ;
d4(1)=d3(1)-d3(2)-6 ; % d4 is X4 d4(2)=d3(1)+d3(2)+6 ; %compute X(t+h) for i=1:n
X(i)=X(i)+h*(d1(i)+h*(d2(i)+h*(d3(i)+h*d4(i)/4)/3)/2); rs(i,k)=X(i) ; end %for i t=t+h ;
chapter6 page 1
end % for k
例題 1: 解聯立常微分方程式
x'(t)=x(t)-y(t) +2 t -t2 –t3 ; y'(t)=x(t)+y(t) - 4t2 +t3 初值x(0)=1, y(0) = 0 , 0 ≦t ≦ 1
a=0; b= 1; m= 10;
X0=[1; 0] ;
Rs =taylor4sys(X0,a,b,m);
Rs = [X0 Rs] %增加初值的資料
Rs =
Columns 1 through 7
1.0000 1.1097 1.2371 1.3796 1.5341 1.6969 1.8639 0 0.1093 0.2347 0.3719 0.5169 0.6654 0.8129 Columns 8 through 11
2.0302 2.1905 2.3389 2.4687 0.9543 1.0845 1.1977 1.2874
t2t3
比較上面的結果與真正解 x(t) = ecos(t) + t ; y(t) = esin(t) –t ;
並畫出圖形
t=0:0.1:1;
sol1= exp(t).*cos(t) + t.^2 sol2= exp(t).*sin(t)- t.^3 tt=0:0.05:1;
xt= exp(tt).*cos(tt) + tt.^2 ; yt= exp(tt).*sin(tt)- tt.^3 ;
plot(Rs(1,:),Rs(2,:), 'r*',xt,yt)
sol1 =
Columns 1 through 7
1.0000 1.1096 1.2371 1.3796 1.5341 1.6969 1.8639 Columns 8 through 11
2.0302 2.1905 2.3389 2.4687 sol2 =
Columns 1 through 7
0 0.1093 0.2347 0.3719 0.5169 0.6654 0.8128 Columns 8 through 11
0.9543 1.0845 1.1977 1.2874
chapter6 page 2
1.41.210.80.60.40.2011.522.5 依據上面的數據我們發現利用
taylor4sys得到的結果, 可正確到小數
點第四位, 此外從圖形上顯示結果亦算蠻好的.
2. rk4sys.m -- 利用Runge-Kutta方法
X(t+h) =X(t) +h /6 (K1+2 K2+2 K3+K4) , K1=F(t,X), K2=F(t+h/2, X+h/2 K1) ,
K2=F(t+h/2, X+h/2 K2) , K4=F(t+h, X+h K3) 來估算聯立常微分方程式 X'(t)=F(X,t)的數值解
顯示rk4sys.m的內容 type rk4sys.m
function rs= rk4sys(X0,a,b,m)
%to return the approximation values of X(t) %by RK4 method, X0 is the initial %t is in[a,b] n=length(X0); rs=zeros(n,m) ;
K1=zeros(n,1) ; K2=zeros(n,1) ; K3=zeros(n,1) ; K4=zeros(n,1) ; t=a ;
X=X0 ;X2=zeros(n,1) ;X3=X2; X4=X2 ; h=(b-a)/m ; for k=1:m
%compute K1,K2,K3and K4 K1=fxsys(t,X); for i=1:n
X2(i)=X(i)+1/2*h*K1(i); end
chapter6 page 3
K2=fxsys(t+h/2,X2); for i=1:n
X3(i)=X(i)+1/2*h*K2(i); end
K3=fxsys(t+h/2,X3); for i=1:n
X4(i)=X(i)+h*K3(i); end
K4=fxsys(t+h,X4); %compute X(t+h) for i=1:n
X(i)=X(i)+h/6*(K1(i)+2*K2(i)+2*K3(i)+K4(i)); rs(i,k)=X(i) ; end %for i t=t+h ; end % for k %
function fx=fxsys(t,X)
%compute K1,K2,K3 and K4 of Runge-Kutta method %
fx=zeros(length(X),1);
fx(1)=X(1)-X(2)+t*(2-t*(1+t)); fx(2)=X(1)+X(2)-t^2*(4-t) ;
例題 2: 解聯立常微分方程式(與例題1 相同)
x'(t)=x(t)-y(t) +2 t -t2 –t3 ; y'(t)=x(t)+y(t) - 4t2 +t3 初值x(0)=1, y(0) = 0 , 0 ≦t ≦ 1
a=0; b= 1; m= 10;
X0=[1; 0] ;
Rs =rk4sys(X0,a,b,m);
Rs = [X0 Rs] %增加初值的資料
Rs =
Columns 1 through 7
1.0000 1.1097 1.2371 1.3796 1.5341 1.6969 1.8639 0 0.1093 0.2347 0.3719 0.5169 0.6654 0.8128 Columns 8 through 11
2.0302 2.1906 2.3389 2.4687 0.9543 1.0845 1.1977 1.2874
從上面的數據我們發現利用rk4sys所得到的結果與taylor4sys的結果
在小數點第四位都一樣 .
3. rk4sysnew.m fxsys.m --
解高階常微分方程式 ,我們是利用解一階聯立常微分方程式的方法. 為了利用
rk4sys.m我們將原程式修正, 把計算函數F(t,X)的部分寫成獨立的 m-function 名為fxsys.m, 而新的程式命名為rk4sysnew.m , 將此二者之內容顯示如下:
type rk4sysnew.m
type fxsys.m
chapter6 page 4
function rs= rk4sysnew(X0,a,b,m)
%to return the approximation values of X(t) %by RK4 method, X0 is the initial %t is in[a,b] n=length(X0); rs=zeros(n,m) ;
K1=zeros(n,1) ; K2=zeros(n,1) ; K3=zeros(n,1) ; K4=zeros(n,1) ; t=a ;
X=X0 ;X2=zeros(n,1) ;X3=X2; X4=X2 ; h=(b-a)/m ; for k=1:m
%compute K1,K2,K3and K4 K1=fxsys(t,X); for i=1:n
X2(i)=X(i)+1/2*h*K1(i); end
K2=fxsys(t+h/2,X2); for i=1:n
X3(i)=X(i)+1/2*h*K2(i); end
K3=fxsys(t+h/2,X3); for i=1:n
X4(i)=X(i)+h*K3(i); end
K4=fxsys(t+h,X4); %compute X(t+h) for i=1:n
X(i)=X(i)+h/6*(K1(i)+2*K2(i)+2*K3(i)+K4(i)); rs(i,k)=X(i) ; end %for i t=t+h ; end % for k
function fx=fxsys(t,X)
%compute values of function F(t,X) %
fx=zeros(length(X),1); fx(1)=X(2);
fx(2)=-3* cos(t)^2 +2;
2
例題 3: 解二階常微分方程式 x''(t)= -3 cos(t) +2
初值 x(0)=0, x'(0) = 0 , 0 ≤t ≤ 1
我們把此二階常微分方程式看成一組一階常微分方程式 令變數 x1 = x ; x2 = x' ;
則原微分方程式可看成 x1' = x2 ; x2' = f(t, x1,x2) = -3 cos2(t) +2 x1(0)=0, x2(0)=0 a=0; b= 1; m= 10;
X0=[0; 0] ;
Rs =rk4sysnew(X0,a,b,m);
Rs = [X0 Rs] %增加初值的資料
chapter6 page 5
Rs =
Columns 1 through 7
0 -0.0050 -0.0196 -0.0430 -0.0737 -0.1099 -0.1491 0 -0.0990 -0.1921 -0.2735 -0.3380 -0.3811 -0.3990 Columns 8 through 11
-0.1888 -0.2259 -0.2577 -0.2811 -0.3891 -0.3497 -0.2804 -0.1820
上面Rs結果第一列表示x(t)在t=0, 0.1, 0.2, …,1.0的數值解
+3/4 cos2(t) –3/4
我們畫出其圖形, 並將上面11個值畫上去 .
此二階微分方程式的真實解為 x(t) = 1/4 t
t=0:0.1:1;
sol= 1/4*t.^2+3/4*(cos(t).^2 -1) plot(t,Rs(1,:), 'r*',t, sol)
sol =
Columns 1 through 7
0 -0.0050 -0.0196 -0.0430 -0.0737 -0.1099 -0.1491 Columns 8 through 11
-0.1888 -0.2259 -0.2577 -0.2811 02
-0.05-0.1-0.15-0.2-0.25-0.3-0.3500.10.20.30.40.50.60.70.80.91 此圖形顯示數值解相當不錯 .
同樣地, 解高階聯立常微分方程式 ,我們也是利用解一階聯立常微分方程式
的方法.
例題4: 解x''= x-y –(3x')2 + (y')3 + 6y'' +2t , y'''=y''-x'+ex –t , 1≤ t ≤1.5 初值x(1)=2, x'(1)=-4, y(1)=-2, y'(1)=7, y''(1)=6
我們令變數 x1 = x, x2=x' , x3= y, x4=y', and x5=y'' 則原方程組可看成 :
chapter6 page 6
x1' = x2 ; x2' = x1- x3 – 9 x22 + x43 +6 x5 +2 t x3' = x4 ; x4' = x5 ; x5' = x5- x2 +ex1 –t
其餘留待給同學當作習題
4. adamsmsys.m – 利用Adams-Moulton (Predictor-Corrector)方法解聯立微分方程組 顯示 adamsmsys.m之內容 type adamsmsys.m
function rs= adamsmsys(X0,X1,X2,X3,a,b,m)
%to return the approximation values of X(t),t is in[a,b] %by Adams-Moulton method, X0,X1,X2,and X3 are the initial úta
n=length(X0); rs=zeros(n,m) ;
rs(:,1)=X1; rs(:,2)=X2; rs(:,3)=X3; K0=zeros(n,1) ; K1=K0 ; K2=K0 ; K3=K0 ; K=K0 ; h=(b-a)/m ; t=a+3*h ; for j=4:m
%compute F(X0),F(X1), F(X2) and F(X3) K0=fxsys(t-3*h,X0); K1=fxsys(t-2*h,X1); K2=fxsys(t-h,X2); K3=fxsys(t,X3); %compute predictor for i=1:n
X(i)=X3(i)+h*(55*K3(i)-59*K2(i)+37*K1(i)-9*K0(i))/24; end
%compute corrector K=fxsys(t+h,X) ; for i=1:n
X(i)=X3(i)+h*(9*K(i)+19*K3(i)-5*K2(i)+K1(i))/24 ; rs(i,j)=X(i) ; end %for i t=t+h ;
X0=X1; X1=X2; X2=X3; X3=X ; end % for j
同樣例題4我們利用此方法觀察執行之結果
例題5: 解x''= x-y –(3x')2 + (y')3 + 6y'' +2t , y'''=y''-x'+ex –t , 1≤ t ≤1.5 初值x(1)=2, x'(1)=-4, y(1)=-2, y'(1)=7, y''(1)=6
我們令變數 x1 = x, x2=x' , x3= y, x4=y', and x5=y'' 則原方程組可看成 :
x1' = x2 ; x2' = x1- x3 – 9 x22 + x43 +6 x5 +2 t x3' = x4 ; x4' = x5 ; x5' = x5- x2 +ex1 –t 所以重新設定函數fxsys.m的內容 type fxsys.m
chapter6 page 7
function fx=fxsys(t,X)
%compute values of function F(t,X) %
fx=zeros(length(X),1); fx(1)=X(2);
fx(2)=X(1)-X(3)-9*X(2)^2+X(4)^3+6*X(5)+2*t; fx(3)=X(4) ; fx(4)=X(5) ;
fx(5)=X(5)-X(2)+exp(X(1))-t ;
因為函數F(t,X)的數值可能很大(小), 為得到較好的解, 我們將間距
h取小,
即在[a, b]區間中取較多的t點(例如200點) .
a=1; b = 1.5 ; m=200 ;
X0= [2 -4 -2 7 6]' ;
%利用rk4sysnew計算X1, X2, and X3 Rs=rk4sysnew(X0,a,1.0075,3);
X1=Rs(:,1); X2=Rs(:,2); X3=Rs(:,3); Rs=adamsmsys(X0,X1,X2,X3,a,b,m);
下面我們取出21組 (x(t),y(t)) 的資料與Matlab 函數
ode45執行的結果比較.
Matlab 提供的ode45函數是利用Runge-Kutta order 4 and 5取不固定間距h 的方法 .
x=zeros(1,20) ; t=x ; y=x ;
for j= 1:20
t(j) = 1+10*j*0.0025 ; x(j) = Rs(1, 10*j) ; y(j) = Rs(3, 10*j) ;
end ;
t=[a t] ; x=[2 x] ;
y=[-2 y]; %增加初值之資料 %
tt = a ;
options = odeset('RelTol', 1e-5, 'AbsTol', 1e-6) ;
[tt , Rst] = ode45('fxsys', [a,b] , X0, options); xx= Rst(:,1); yy=Rst(:,3) ;
plot(tt,xx,tt,yy,t,x,'r*',t,y,'g*') text(1.2,2.4,'x(t)','FontSize',16)
text(1.2, -1.2, 'y(t)','FontSize',16)
chapter6 page 8
876543210-1-211.051.11.15x(t)y(t)1.21.251.31.351.41.451.5 圖形顯示我們利用
adamsmsys函數執行的結果還算不錯.
剛才提到Matlab 提供的ode45函數是取不固定間距h的方法, 同樣地,
在Adams-Moulton (predictor-corrector) 的方法中, 若考慮局部誤差
τi = (19/270 ) |Xi –Xi(p)| / |Xi| 小於設定的ε則維持原來間距h否則間距h變 成h/2 , 再重新計算相關的點.
chapter6 page 9