一阶微分方程式初值问题.docx
《一阶微分方程式初值问题.docx》由会员分享，可在线阅读，更多相关《一阶微分方程式初值问题.docx（14页珍藏版）》请在冰点文库上搜索。
一阶微分方程式初值问题
第六章一階微分方程式初值問題
所謂一階微分方程式初值問題,即求滿足下列微分方程式的解
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
升级会员 