一阶微分方程式初值问题.docx
《一阶微分方程式初值问题.docx》由会员分享,可在线阅读,更多相关《一阶微分方程式初值问题.docx(12页珍藏版)》请在冰豆网上搜索。
一阶微分方程式初值问题
第六章一階微分方程式初值問題
所謂一階微分方程式初值問題,即求滿足下列微分方程式的解
x'(t)=f(t,x),x(a)=x0
一般分為單一點(singlestep)計算與多點(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:
如果你有安裝MatlabNotebook要執行下列inputcells(綠色指令敘述)
之前必須先執行上面的cell–[path(…)]
藍色的內容是Matlab[outputcells]
註2:
如果m-file之內有定義函數,請記得改寫你要執行的.
1.taylor4.m–
顯示taylor4.m的內容
typetaylor4.m
functiony=taylor4(x0,a,b,m)
%toreturntheapproximationvaluesofx(t)
%byTaylorseriesmethodoforder4,x0istheinitial
%tisin[a,b]
y=[];
n=1;
t=a;
x=x0;
h=(b-a)/m;
%fprintf('taylor4\n')
%fprintf('ntx\n')
fork=1:
m
%Forgivenx'(t)=f(t,x)=1+x^2+t^3
%computederivativesofx',x'',x'''andx4att
d1=1+x^2+t^3;%d1isx'
d2=2*x*d1+3*t^2;%d2isx''
d3=2*x*d2+2*d1^2+6*t;%d3isx'''
d4=2*x*d3+6*d1*d2+6;%d4isx""
%computex(t+h)
x=x+h*(d1+h*(d2+h*(d3+h*d4/4)/3)/2);
t=t+h;
%fprintf('%3d%10.6f%10.6f\n',k,t,x)
y=[y;ktx];
end%fork
2.rk4.m–
顯示rk4.m的內容
typerk4.m
functionrs=rk4(x0,a,b,m)
%toreturntheapproximationvaluesofx(t)
%byRK4method,x0istheinitial
%tisin[a,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('rk4\n')
%fprintf('ntx\n')
fork=1:
m
%computeK1,K2,K3andK4
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);
x4=x+h*K3;
K4=fx(t+h,x4);
%computex(t+h)
x=x+h*(K1+2*K2+2*K3+K4)/6;
t=t+h;
%fprintf('%3d%10.6f%10.6f\n',k,t,x)
rs=[rs;ktx];
end%fork
%
functiony=fx(t,x)
%computevaluesoffunctionf(t,x)
%
y=1+x^2+t^3;
例題1:
比較Taylor與RK的方法解微分方程式
x'(t)=1+x2+t3;x(0)=-4;1≦t≦2
ty4rk4
taylorvsrk4
nttaylorrk4
11.050000-3.246802-3.245493
21.100000-2.696215-2.694804
31.150000-2.268329-2.267053
41.200000-1.918709-1.917595
51.250000-1.620467-1.619494
61.300000-1.356111-1.355251
71.350000-1.113464-1.112692
81.400000-0.883455-0.882749
91.450000-0.658806-0.658147
101.500000-0.433177-0.432548
111.550000-0.200533-0.199920
121.6000000.0454230.046037
131.6500000.3118710.312505
141.7000000.6076840.608360
151.7500000.9446450.945400
161.8000001.3394921.340382
171.8500001.8176151.818749
181.9000002.4204022.422010
191.9500003.2212943.223954
202.0000004.3657344.371224
Multi-stepmethod
例題2:
比較Adams-Bashforth與Milne的方法解微分方程式
x'(t)=-6x+6;x(0)=2;0≦t≦1
3.Adams-Bashforth--
顯示adbh.m的內容
typeadbh.m
functionrs=adbh(x0,a,b,m)
%toreturntheapproximationvaluesofx(t)
%by4th-orderAdams-Bashforthmethod,x0istheinitial
%tisin[a,b]
rs=[];
x=zeros(m,1);t=x;
h=(b-a)/m;
inl=[00x0];
%obtaintheinitialsbyrk4
rs=rk4adbh(x0,a,a+3*h,3);
rs=[inl;rs];
k=rs(:
1);t=rs(:
2);x=rs(:
3);
%fprintf('rk4\n')
%fprintf('ntx\n')
%computex(t+h)
fork=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;
t(k+1)=t(k)+h;
end
ext=inline('1+exp(-6*t)','t');
xt=feval(ext,t);
%fork=1:
m+1
%fprintf('%3d%10.6f%10.6f\n',k,t(k),x(k))
%end
k=(1:
m+1)';
rs=[ktxxt];
%
functiony=f(x)
%computevaluesoffunctionf(t,x)
%
y=-6.*x+6;
4.milne–
顯示milne.m的內容
typemilne.m
functionrs=milne(x0,a,b,m)
%toreturntheapproximationvaluesofx(t)
%by4th-orderMilne'smethod,x0istheinitial
%tisin[a,b]
rs=[];x=zeros(m,1);t=x;
h=(b-a)/m;
inl=[00x0];
%obtaintheinitialsbyrk4usingthesamef(t,x)as
%Adams-Bashforth
rs=rk4adbh(x0,a,a+3*h,3);
rs=[inl;rs];
k=rs(:
1);t=rs(:
2);x=rs(:
3);
%fprintf('rk4\n')
%fprintf('ntx\n')
%computex(t+h)
fork=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);
%fork=1:
m+1
%fprintf('%3d%10.6f%10.6f\n',k,t(k),x(k))
%end
k=(1:
m+1)';
rs=[ktxxt];
%
functiony=f(x)
%computevaluesoffunctionf(t,x)
%
y=-6.*x+6;
adbhmilne
Adams-BashforthvsMilne
ntAdams-BFMilneexactsolution
10.0000002.0000002.0000002.000000
20.1000001.5494001.5494001.548812
30.2000001.3018401.3018401.301194
40.3000001.1658311.1658311.165299
50.4000001.0998331.0971031.090718
60.5000001.0515761.0437561.049787
70.6000001.0424331.0441831.027324
80.7000001.0051290.9747801.014996
90.8000001.0354191.1027911.008230
100.9000000.9666380.7884221.004517
111.0000001.0695571.5052921.002479
從上面數據顯示Adam-bashforth穩定度較好
5.trapzoidnw.m–
顯示trapzoidnw.m的內容
typetrapzoidnw.m
function[msg,rs]=trapzoidnw(x0,a,b,m,tor)
%toapproximatethesolutionofx'=f(t,x)
%tisin[a,b],x(a)=x0byTrapezoidalwith
%NewtonIteration,iterationlimitis20
rs=[];
x=zeros(m,1);t=x;
h=(b-a)/m;
x
(1)=x0;t
(1)=a;
%computex(t+h)
fork=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);
iferr<=tor%acceptw
t(k+1)=th;x(k+1)=w;
else
j=j+1;
w0=w;
end;
end;%while
if(j>20)
msg='Newtoniterationfails';
break;
else
msg='Newtoniterationsuccess';
end;%if
end;%fork
ext=inline('t-exp(-5*t)','t');
xt=feval(ext,t);
%fork=1:
m+1
%fprintf('%3d%10.6f%10.6f\n',k,t(k),x(k))
%end
rk=(1:
m+1)';
rs=[rktxxt];
functiony=f(t,x)
%computevaluesoffunctionf(t,x)
%
y=5*exp(5*t)*(x-t)^2+1;
functiony=fy(t,x)
%computevaluesoffunctionDf(t,x)/Dx
%
y=10*exp(5*t)*(x-t);
例題3:
利用隱性梯形法(ImplicitTrapezoidalmethod)解微分方程式
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('----TrapezoidalwithNewtoniterationvsExactsolution----\n\n')
fprintf('ntTrapezoidalexactsolution\n')
fork=1:
m+1
fprintf('%3d%10.4f%10.6f%10.6f\n',k,rs(k,2),rs(k,3),rs(k,4))
end
----TrapezoidalwithNewtoniterationvsExactsolution----
ntTrapezoidalexactsolution
10.0000-1.000000-1.000000
20.1000-0.501080-0.506531
30.2000-0.162742-0.167879
40.30000.0806070.076870
50.40000.2671430.264665
60.50000.4194900.417915
70.60000.5511930.550213
80.70000.6704050.669803
90.80000.7820530.781684
100.90000.8891160.888891
111.00000.9933990.993262
由於採用牛頓法須要計算函數的微分較為複雜,
故我們利用Eulermethod當做predcitor,Trapezodialmethod
當做corrector,得到下列結果.
a=0;b=1;x0=-1;m=10;
rs=trapzoidpc(x0,a,b,m);
fprintf('----Trapezoidal(P-C)vsExactsolution----\n\n')
fprintf('ntTrapezoidal(P-C)exactsolution\n')
fork=1:
m+1
fprintf('%3d%10.4f%10.6f%10.6f\n',k,rs(k,2),rs(k,3),rs(k,4))
end
----Trapezoidal(P-C)vsExactsolution----
ntTrapezoidal(P-C)exactsolution
10.0000-1.000000-1.000000
20.1000-0.546955-0.506531
30.2000-0.212491-0.167879
40.30000.0399390.076870
50.40000.2374650.264665
60.50000.3991070.417915
70.60000.5377030.550213
80.70000.6616940.669803
90.80000.7765210.781684
100.90000.8856450.888891
111.00000.9912400.993262
其誤差比上述之結果稍大.同學請比較其他的方法.
6.eulerdynamic.m--
雖然利用Euler的方法:
wj+1=wj+h*f(t,wj),j>=0,w0=x0
估算微分方程式的解x'(t)=f(t,x),x(a)=x0精確度不佳,但是我們利用此方法
所產生之數據,來觀察動態系統(dynamicsystem)在平衡狀態可能發生之現象:
perioddoublingandchaos.
例題:
微分方程式
p'(t)=10p(1-p);p(0)=0.1
顯示eulerdynamic.m的內容
typeeulerdynamic.m
%toplot30approximationvaluesofx(t)after200iterationsbyEuler's
%methodtodemothephenomenaofperioddoubleingandchaos
%ofdynamicalsystem
y=[];
m=230;
pos=[1020360300];
figure('Position',pos)
holdon
forh=0.18:
0.001:
0.3
t=0;x=0.1;
fork=1:
m
%Forgivenx'(t)=f(t,x)=10x(1-x)
%computex(t+h)
x=(1+10*h)*x-(10*h)*x^2;
t=t+h;
ifk>200
y=[y;ktx];
ifrem(k,4)==0
plot(h,x,'y.')
elseifrem(k,4)==1
plot(h,x,'r.')
elseifrem(k,4)==2
plot(h,x,'g.')
else
plot(h,x,'b.')
end
end
end%fork
end%h
holdoff
執行
eulerdynamic