1、第六章 一階微分方程式初值問題 所謂一階微分方程式初值問題, 即求滿足下列微分方程式的解 x'(t) = f(t, x), x(a) = x0 一般分為單一點(single step)計算與多點 (multi-step)計算的方法. 本章的m-file 如下: 1. taylor4.m 2. rk4.m 3. adbh.m (Adams-Bashforth) 4. milne.m 5. trapzoidnw.m 6. eulerdynamic.m 將須要的m-file之檔案夾加入搜尋路徑中 path('d:numerical', path) 註1: 如果你有安裝Matl

2、ab Notebook 要執行下列input cells (綠色指令敘述) 之前必須先執行上面的cell path () 藍色的內容是Matlab output cells 註2: 如果 m-file之內有定義函數, 請記得改寫你要執行的 . 1. taylor4.m 顯示taylor4.m的內容 type taylor4.m function y=taylor4(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

3、,b y=; n=1; t=a ; x=x0 ; h=(b-a)/m ; % fprintf('taylor4n') % fprintf(' n t xn') for k=1:m%For given x'(t)= f(t,x)=1+x2+t3%compute derivatives of x', x'', x''' and x4 at t d1=1+x2+t3 ; % d1 is x' d2=2*x*d1 + 3*t2; % d2 is x'' d3=2*x*d2 + 2*d12 +

4、 6*t ; % d3 is x''' d4=2*x*d3 + 6*d1*d2+ 6 ; % d4 is x"" %compute x(t+h) x=x+h*(d1+h*(d2+h*(d3+h*d4/4)/3)/2); t=t+h ; % fprintf('%3d %10.6f %10.6fn',k ,t ,x) y=y; k t x; end % for k 2. rk4.m 顯示rk4.m的內容 type rk4.m function rs= rk4(x0,a,b,m)%to return the approximation va

5、lues of x(t)%by RK4 method, x0 is the initial%t is ina,b rs= ; K1=0 ; K2=0 ; K3=0 ; K4=0; t=a ; x=x0 ;x2=0 ;x3=x2; x4=x2 ; h=(b-a)/m ; % fprintf('rk4n') % fprintf(' n t xn') for k=1:m%compute K1,K2,K3and K4 K1=fx(t,x); x2=x+1/2*h*K1; K2=fx(t+h/2,x2); x3=x+1/2*h*K2; K3=fx(t+h/2,x3); x

6、4=x+h*K3; K4=fx(t+h,x4);%compute x(t+h) x=x+h*(K1+2*K2+2*K3+K4)/6; t=t+h ; % fprintf('%3d %10.6f %10.6fn',k ,t ,x) rs=rs; k t x; end % for k %function y=fx(t,x)%compute values of function f(t,x)% y= 1+x2+t3; 例題 1: 比較Taylor與 RK的方法解微分方程式 x'(t)=1+x2 +t3 ; x(0)= -4 ; 1 t 2 ty4rk4 taylor vs r

7、k4n t taylor rk4 1 1.050000 -3.246802 -3.245493 2 1.100000 -2.696215 -2.694804 3 1.150000 -2.268329 -2.267053 4 1.200000 -1.918709 -1.917595 5 1.250000 -1.620467 -1.619494 6 1.300000 -1.356111 -1.355251 7 1.350000 -1.113464 -1.112692 8 1.400000 -0.883455 -0.882749 9 1.450000 -0.658806 -0.658147 10 1

8、.500000 -0.433177 -0.432548 11 1.550000 -0.200533 -0.199920 12 1.600000 0.045423 0.046037 13 1.650000 0.311871 0.312505 14 1.700000 0.607684 0.608360 15 1.750000 0.944645 0.945400 16 1.800000 1.339492 1.340382 17 1.850000 1.817615 1.818749 18 1.900000 2.420402 2.422010 19 1.950000 3.221294 3.223954

9、20 2.000000 4.365734 4.371224 Multi-step method 例題 2: 比較Adams-Bashforth與 Milne的方法解微分方程式 x'(t)= -6x+ 6 ; x(0)= 2 ; 0 t 1 3. Adams-Bashforth - 顯示adbh.m的內容 type adbh.m function rs= adbh(x0,a,b,m)%to return the approximation values of x(t)%by 4th-order Adams-Bashforth method, x0 is the initial%t is

10、ina,b rs= ; x=zeros(m,1); t=x; h=(b-a)/m ; inl=0 0 x0; %obtain the initials by rk4 rs= rk4adbh(x0,a,a+3*h, 3); rs=inl; rs ; k=rs(:,1); t=rs(:,2); x=rs(:,3); % fprintf('rk4n') % fprintf(' n t xn')%compute x(t+h) for k=4:m x(k+1)= x(k)+h*(55*f(x(k)-59*f(x(k-1)+37*f(x(k-2)-9*f(x(k-3)/24

11、 ; t(k+1)=t(k)+h ; end ext =inline(' 1+exp(-6*t)','t'); xt=feval(ext,t) ;% for k=1:m+1% fprintf('%3d %10.6f %10.6fn',k ,t(k) ,x(k)% end k=(1:m+1)' rs=k t x xt; %function y=f(x)%compute values of function f(t,x)% y= -6.*x + 6; 4. milne 顯示milne.m的內容 type milne.m function rs

12、= milne(x0,a,b,m)%to return the approximation values of x(t)%by 4th-order Milne's method, x0 is the initial%t is ina,b rs= ; x=zeros(m,1); t=x; h=(b-a)/m ; inl=0 0 x0; %obtain the initials by rk4 using the same f(t,x) as %Adams-Bashforth rs= rk4adbh(x0,a,a+3*h, 3); rs=inl; rs ; k=rs(:,1); t=rs(:

13、,2); x=rs(:,3); % fprintf('rk4n') % fprintf(' n t xn')%compute x(t+h) for k=4:m x(k+1)= x(k-3)+4*h*(2*f(x(k)-f(x(k-1)+2*f(x(k-2)/3 ; t(k+1)=t(k)+h ; end ext =inline(' 1+exp(-6*t)','t'); xt=feval(ext,t) ;% for k=1:m+1% fprintf('%3d %10.6f %10.6fn',k ,t(k) ,x(k)

14、% end k=(1:m+1)' rs=k t x xt; %function y=f(x)%compute values of function f(t,x)% y= -6.*x + 6; adbhmilne Adams-Bashforth vs Milnen t Adams-BF Milne exact solution 1 0.000000 2.000000 2.000000 2.000000 2 0.100000 1.549400 1.549400 1.548812 3 0.200000 1.301840 1.301840 1.301194 4 0.300000 1.16583

15、1 1.165831 1.165299 5 0.400000 1.099833 1.097103 1.090718 6 0.500000 1.051576 1.043756 1.049787 7 0.600000 1.042433 1.044183 1.027324 8 0.700000 1.005129 0.974780 1.014996 9 0.800000 1.035419 1.102791 1.008230 10 0.900000 0.966638 0.788422 1.004517 11 1.000000 1.069557 1.505292 1.002479 從上面數據顯示Adam-

16、bashforth穩定度較好5. trapzoidnw.m 顯示trapzoidnw.m的內容 type trapzoidnw.m function msg, rs= trapzoidnw(x0, a,b, m, tor)%to approximate the solution of x'=f(t,x)%t is in a,b, x(a)=x0 by Trapezoidal with%Newton Iteration , iteration limit is 20 rs= ; x=zeros(m,1); t=x; h=(b-a)/m ; x(1)=x0; t(1)=a; %comput

17、e x(t+h) for k=1:m c=x(k)+h*f(t(k),x(k)/2 ; w0=x(k) ; j=1; err = 1.0 ; th=t(k)+h ; while(j<=20 & err > tor) w=w0 -(w0-h/2*f(th,w0)-c)/(1-h/2*fy(th,w0) ; err=abs(w- w0) ; if err<= tor %accept w t(k+1)=th; x(k+1)=w ; else j=j+1 ; w0=w ; end ; end; %while if (j>20) msg='Newton itera

18、tion fails' break ; else msg='Newton iteration success' ; end ; %if end ; % for k ext =inline('t-exp(-5*t)','t'); xt=feval(ext,t) ;% for k=1:m+1% fprintf('%3d %10.6f %10.6fn',k ,t(k) ,x(k)% end rk=(1:m+1)' ; rs=rk t x xt; function y=f(t,x)%compute values of fu

19、nction f(t,x)% y= 5*exp(5*t)*(x-t)2+1;function y=fy(t,x)%compute values of function Df(t,x)/Dx% y= 10*exp(5*t)*(x-t); 例題 3: 利用隱性梯形法(Implicit Trapezoidal method)解微分方程式 x'(t)= 5e5t(x-t)2 +1 ; x(0)= -1 ; 0 t 1 a=0; b= 1; x0= -1; m=10; tor = 10-6 ; msg, rs = trapzoidnw(x0, a, b, m, tor) ; fprintf(&#

20、39;-Trapezoidal with Newton iteration vs Exact solution-nn') fprintf(' n t Trapezoidal exact solutionn') for k=1:m+1 fprintf('%3d %10.4f %10.6f %10.6fn',k,rs(k,2), rs(k,3), rs(k,4) end -Trapezoidal with Newton iteration vs Exact solution- n t Trapezoidal exact solution 1 0.0000 -

21、1.000000 -1.000000 2 0.1000 -0.501080 -0.506531 3 0.2000 -0.162742 -0.167879 4 0.3000 0.080607 0.076870 5 0.4000 0.267143 0.264665 6 0.5000 0.419490 0.417915 7 0.6000 0.551193 0.550213 8 0.7000 0.670405 0.669803 9 0.8000 0.782053 0.781684 10 0.9000 0.889116 0.888891 11 1.0000 0.993399 0.993262 由於採用牛

22、頓法須要計算函數的微分較為複雜, 故我們利用Euler method 當做predcitor, Trapezodial method 當做corrector, 得到下列結果 . a=0; b= 1; x0= -1; m=10; rs = trapzoidpc(x0, a, b, m) ; fprintf('-Trapezoidal(P-C) vs Exact solution-nn') fprintf(' n t Trapezoidal(P-C) exact solutionn') for k=1:m+1 fprintf('%3d %10.4f %10.

23、6f %10.6fn',k,rs(k,2), rs(k,3), rs(k,4) end -Trapezoidal(P-C) vs Exact solution- n t Trapezoidal(P-C) exact solution 1 0.0000 -1.000000 -1.000000 2 0.1000 -0.546955 -0.506531 3 0.2000 -0.212491 -0.167879 4 0.3000 0.039939 0.076870 5 0.4000 0.237465 0.264665 6 0.5000 0.399107 0.417915 7 0.6000 0.537703 0.550213 8 0.7000 0.661694 0.669803 9 0.8000 0.776521 0.781684 10 0.9000 0.885645 0.888891 11 1.0000 0.991240 0.993262 其誤差比上述之結果稍大. 同學請比較其他的方法 . 6. eulerdynamic.m - 雖然利用Euler的方法: wj+1 = wj+ h*f(t,wj) , j>=0, w0 = x0估算微分方程式的解 x'(t) = f