Matlab-求解化工常微分方程和偏微分方程課件

Matlab求解化工常微分方程和偏微分方程方利國Matlab求解化工常微分方程和偏微分方程方利國1Matlab求解化工常微分方程和偏微分方程1、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命令1.2初值問題求解1.3邊值問題求解1.4加權(quán)問題求解（自學內(nèi)容）2、偏微分方程（組）求解2.1問題描述及一維動態(tài)PDE方程求解2

.2二維求解Matlab求解化工常微分方程和偏微分方程1、常微分方程（21、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命令

常微分方程：（初值問題）常微分方程：（兩點邊值問題）

1、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命31、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命令

微分方程組：

1、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命41、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命令

高級微分方程：

1、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命51、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命令Matlab調(diào)用命令：ODE45:4-5階龍格庫塔法（非剛性）ODE23:2-3階龍格庫塔法（非剛性）ODE113:可變D-B-M法（非剛性）ODE15S:基于數(shù)值差分的可變階方法法（剛性）ODE23S、ODE23t、ODE23tb（剛性）

1、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命61、常微分方程（組）求解通用調(diào)用格式：[x,y]=ode***(@odefun,xspan,y0,<option,p1,p2..>)X:自變量向量，在實際調(diào)用時取名不一定要用x，也可以用其他名稱，只要前后一致即可。Y:應(yīng)變量向量，在實際調(diào)用時取名不一定要用Y，也可以用其他名稱，只要前后一致即可。***:根據(jù)不同的問題調(diào)用不同格式，如45，23s@odefun:自定義函數(shù)的函數(shù)名，該函數(shù)為Xspan：自變量的積分限，[xa,xb],也可以是離散點，[x0,x1,x2,…xf]y0：應(yīng)變量向量的初值<>:可以沒有該選項，如有，具體應(yīng)用見下面的實際例子

1、常微分方程（組）求解通用調(diào)用格式：71.2初值問題求解例1：已知某高溫物體其溫降過程符合以下規(guī)律，其中溫度T的單位為K，時間

的單位為分鐘，零時刻高溫物體的溫度為2000K，以1分鐘作為時間步長，請計算零時刻以后每隔1分鐘至170分鐘的溫度。單個微分方程1.2初值問題求解例1：已知某高溫物體其溫降過程符合以下8functionxODEs%鐵球從2000K降溫曲線,在7.0版本上調(diào)試通過%由華南理工大學方利國編寫，2012年2月29日%歡迎讀者調(diào)用，如有問題請告知lgfang@clearall；clcy0=[2000];[x1,y1]=ode45(@f,[0:1:170],y0);%0到170分鐘，每分鐘一個計算點[x2,y2]=ode23(@f,[0:1:170],y0);plot(x1,y1,'r-')xlabel('時間，M')ylabel('溫度，K')holdondisp('Resultsbyusingode45()：')disp('xy(1)')disp([x1y1])disp('Resultsbyusingode23()：')disp('xy(2)')disp([x2y2])plot(x2,y2,'k:')%------------------------------------------------------------------functiondy=f(x,y)%定義降溫速率的微分方程%dy=0.04*y(1)-100;dy=-0.04*exp(0.001*(y(1)-300))*(-300+y(1));functionxODEs9Resultsbyusingode45()：

xy(1)1.0e+003*Columns1through1300.01000.02000.03000.04000.05000.06000.07000.08000.09000.10000.11000.12002.00000.87880.61330.48800.41810.37620.34980.33280.32180.31450.30970.30650.3043Columns14through180.13000.14000.15000.16000.17000.30290.30190.30130.30090.3006Resultsbyusingode23()：

xy(2)1.0e+003*Columns1through1300.01000.02000.03000.04000.05000.06000.07000.08000.09000.10000.11000.12002.00000.87790.61240.48730.41760.37560.34930.33230.32130.31400.30920.30610.3041Columns14through180.13000.14000.15000.16000.17000.30270.30180.30120.30080.3005Resultsbyusingode45()：10Matlab-求解化工常微分方程和偏微分方程課件111.2初值問題求解該問題相當與一個自變量，兩個應(yīng)變量問題，已知初值及微分表達式，可以利用ODE45求解。微分方程組求解1.2初值問題求解該問題相當與一個自變量，兩個應(yīng)變量問題12程序代碼functionuvDEs%微生物消亡問題計算,在7.0版本上調(diào)試通過%由華南理工大學方利國編寫，2012年3月12日%歡迎讀者調(diào)用，如有問題請告知lgfang@clearall;Clcy0=[1.61.2];[x1,y1]=ode45(@f,[0:0.1:10],y0);%0到3分鐘，每0.1分鐘一個計算點u=y1(:,1);v=y1(:,2);plot(x1,u,'r-')xlabel('時間，M')ylabel('微生物濃度')holdonplot(x1,v,'k:')disp('Resultsbyusingode45()：')disp('xuv')disp([x1y1])%------------------------------------------------------------------functiondy=f(x,y)%定義降溫速率的微分方程f1=0.09*y(1)*(1-y(1)/20)-0.45*y(1)*y(2);f2=0.06*y(2)*(1-y(2)/15)-0.001*y(1)*y(2);dy=[f1;f2];程序代碼functionuvDEs131.2初值問題求解

例3：當X較大時，兩種方法計算結(jié)果有較大不同，為什么？單個微分方程有零點問題？1.2初值問題求解例3：當X較大時，兩種方法計算結(jié)果有14functionL43ODEs%在7.0版本上調(diào)試通過%由華南理工大學方利國編寫，2012年2月29日%歡迎讀者調(diào)用，如有問題請告知lgfang@clearallclcy0=[1];[x1,y1]=ode45(@f,[0:0.05:10],y0);%0到10，每0.05間隔一個計算點[x2,y2]=ode23(@f,[0:0.05:10],y0);%0到10，每0.05間隔一個計算點plot(x1,y1,'r-')xlabel('x')ylabel('y')holdondisp('Resultsbyusingode45()：')disp('xy(1)')disp([x1y1])disp('Resultsbyusingode23()：')disp('xy(2)')disp([x2y2])plot(x2,y2,'b--')%------------------------------------------------------------------functiondy=f(x,y)%定義微分方程dy=y^2*cos(x);functionL43ODEs15計算值xy(1)Columns1through1300.05000.10000.15000.20000.25000.30000.35000.40000.45000.50000.55000.60001.00001.05261.11091.17571.24791.32871.41951.52181.63781.76981.92102.09512.2970Columns14through260.65000.70000.75000.80000.85000.90000.95001.00001.05001.10001.15001.20001.25002.53292.81073.14113.53814.02064.61525.35986.30787.54359.192811.461414.728319.5991Columns27through291.30001.35001.400027.428341.203068.6630計算值xy(1)16高階微分方程求解求解思路：將高階微分方程通過變量轉(zhuǎn)換，轉(zhuǎn)變成一級微分方程組進行求解。例4：高階微分方程求解求解思路：將高階微分方程通過變量轉(zhuǎn)換，轉(zhuǎn)變成17高階微分方程求解程序及解Columns1through1300.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.10001.200000.10020.20150.30520.41290.52630.64760.77890.92301.08271.26141.46271.69081.00001.00531.02271.05441.10261.16981.25881.37231.51361.68611.89352.13982.4296Columns14through211.30001.40001.50001.60001.70001.80001.90002.00001.95022.24612.58412.97053.41233.91714.49365.15132.76773.15953.61104.12864.71965.39186.15417.0159高階微分方程求解程序及解Columns1through18剛性方程求解有些微分組的系數(shù)變化很大，這時用ODE45就很難收斂求解，這時可用專門解決此類微分方程的ODE23S來求解，需要注意的是在解的圖像繪制時，也需要考慮數(shù)值的波動幅度很大，需要引入對數(shù)坐標。例5：

dy1/dx=-0.03*y1+1e4*y2*y3dy2/dx=0.03*y1-2e4*y2*y3-5e7*y2^2dy3/dx=5e7*y2^2y1(0)=1,y2(0)=0,y3(0)=0剛性方程求解有些微分組的系數(shù)變化很大，這時用ODE45就很難19functiongangxinDEs%剛性問題計算,在7.0版本上調(diào)試通過%由華南理工大學方利國編寫，2012年3月13日%歡迎讀者調(diào)用，如有問題請告知lgfang@%dy1=-0.03*y1+1e4*y2*y3%dy2=0.03*y1-2e4*y2*y3-5e7*y2^2%dy3=5e7*y2^2clearallclcfigurexspan=[06*logspace(-6,6)];y0=[100];[x1,y1]=ode15s(@f,xspan,y0);%用ode45計算剛性方程，可能有問題u=y1(:,1);v=1e4*y1(:,2);w=y1(:,3);semilogx(x1,u,'r-','linewidth',2)xlabel('x')ylabel('1e4*v')holdonsemilogx(x1,v,'k:','linewidth',2)holdonsemilogx(x1,w,'g-','linewidth',2)gridaxis([10^(-10)10^10-0.21.2])legend('u','v','w')disp('Resultsbyusingode45()：')disp('xuvw')formatlongdisp([x1y1])%------------------------------------------------------------------functiondy=f(x,y)%定義微分方程f1=-0.03*y(1)+1e4*y(2)*y(3);f2=0.03*y(1)-2e4*y(2)*y(3)-5e7*y(2)^2;f3=5e7*y(2)^2;dy=[f1;f2;f3];functiongangxinDEsdisp('x201.3邊值問題求解邊值問題相對于初值問題而言，多了一個端點的約束，如果在高階或微分方程組中端點約束過多，微分方程組可能無解，端點約束有一定限制?？梢酝ㄟ^建立離散的方程組，再利用ODE45進行求解，但可以利用MATLAB的專用工具求解最好。下面介紹ODE-BVPs的求解器，主要有bvpinit,bvp4c,deval,solinit等。通過實際例子介紹這些內(nèi)部函數(shù)的功能。1.3邊值問題求解邊值問題相對于初值問題而言，多了一個端點211.3邊值問題求解solinit=bvpinit（x,yinit):產(chǎn)生在初始網(wǎng)格上的初始解，以便bvp4c調(diào)用，其中x為自變量網(wǎng)格，yinit為對應(yīng)函數(shù)的初值。sol=bvp4c(@odefun,@BCfun,solinit,<option,p1,p2..>)@Bcfun:為定義邊界條件方程，Bcfun(ya,yb),其中ya、yb分別表示左右邊界。其他符號意義同上。deval(sol,xint):計算任意點處的函數(shù)值。1.3邊值問題求解solinit=bvpinit（x,yi221.3邊值問題求解

例6：1.3邊值問題求解例6：23程序functionBVP4c1%求解兩點邊值問題的示例在7.0版本上調(diào)試通過%由華南理工大學方利國編寫，2012年3月13日%歡迎讀者調(diào)用，如有問題請告知lgfang@clearallclca=0;b=10;solinit=bvpinit(linspace(a,b,101),[00]);sol=bvp4c(@ODEfun,@BCfun,solinit);formatshortx=[0:0.1:10];y=deval(sol,x);y1=y(1,:)y2=y(2,:)plot(x,y1,'r-')xlabel('x')ylabel('y')holdongridplot(x,y2,'k:')legend('y','dy')disp('Resultsbyusingbvp4c：')disp('xydy')disp([x1y1])%------------------------------------------------------------------functiondy=ODEfun(x,y)f1=y(2);f2=0.05*(1+x^2)*y(1)+2;dy=[f1;f2];%------------------------------------------------------------------functionbc=BCfun(ya,yb)bc=[ya(1)-40;yb(1)-80];Resultsbyusingbvp4c：

xydyColumns1through1300.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.10001.200039.999838.172236.384034.634732.924431.253129.621428.029826.479124.970223.503822.081020.7025-18.4739-18.0779-17.6872-17.2985-16.9088-16.5158-16.1176-15.7125-15.2996-14.8781-14.4476-14.0080-13.5597程序functionBVP4c1disp('Results24計算結(jié)果計算結(jié)果251.3邊值問題求解

例7：1.3邊值問題求解例7：26主要程序

functiondy=ODEfun(x,y)f1=y(2);f2=4*y(1);dy=[f1;f2];%------------------------------------------------------------------functionbc=BCfun(ya,yb)bc=[ya(1)-1;yb(1)-exp(4)];主要程序

functiondy=ODEfun(x,y)27functionBVP4c2%求解兩點邊值問題的示例在7.0版本上調(diào)試通過%由華南理工大學方利國編寫，2012年3月13日%歡迎讀者調(diào)用，如有問題請告知lgfang@clearallclca=0;b=2solinit=bvpinit(linspace(a,b,21),[00]);sol=bvp4c(@ODEfun,@BCfun,solinit);formatshortx=[0:0.1:2];y=deval(sol,x);y1=y(1,:);y2=y(2,:);y3=exp(2*x);plot(x,y1,'r-','linewidth',2)xlabel('x')ylabel('y')holdongridplot(x,y2,'k:','linewidth',2)holdonplot(x,y3,'b:','linewidth',2)legend('y','dy','real-y')disp('Resultsbyusingbvp4c：')disp('xydy')disp([x;y1;y2])%------------------------------------------------------------------functiondy=ODEfun(x,y)f1=y(2);f2=4*y(1);dy=[f1;f2];%------------------------------------------------------------------functionbc=BCfun(ya,yb)bc=[ya(1)-1;yb(1)-exp(4)];functionBVP4c2282、偏微分方程（組）求解

2.1問題描述2

.2一維動態(tài)PDE方程求解

2

.3二維穩(wěn)態(tài)PDE方程求解(自學）

2、偏微分方程（組）求解2.1問題描述29問題描述當A，B，C為常數(shù)時，稱為擬線性偏微分方程，當A，B，C滿足不同條件時，分為三種不同的類型：不同類型的方程，MATLAB求解時，采用不同或相同的內(nèi)置函數(shù)或工具箱問題描述當A，B，C為常數(shù)時，稱為擬線性偏微分方程，當A，B30PDE求解數(shù)學模型通式m=0，表示平板，m=1表示圓柱，m=2表示球形，f項表示通量項，,s項表示源項，c項為對角陣，元素必須大于等于0才可以求解。PDE求解數(shù)學模型通式m=0，表示平板，m=1表示圓柱，31調(diào)用方法sol=pdepe(m,@pdefun,,@iCfun@BCfun,xspan,tspan,<option,p1,p2..>)@Bcfun:為定義邊界條件方程@iCfun:為定義初始條件方程。其他符號意義同上，注意邊界條件必須寫成以上形式：調(diào)用方法sol=pdepe(m,@pdefun,,@iCf32應(yīng)用策略Pdepe內(nèi)部函數(shù)具體應(yīng)用時，需將實際的偏微分方程對照標準模型，確定c、f、s函數(shù)的具體形式及邊界條件p、q的具體形式及m值。蒸汽入口流體入口，u，t0冷凝液出口流體出口，u，tL應(yīng)用策略Pdepe內(nèi)部函數(shù)具體應(yīng)用時，需將實際的偏微分方程33套管動態(tài)傳熱問題原方程：代入具體數(shù)據(jù)，并對長度歸一化無量綱處理后，得到以下方程：套管動態(tài)傳熱問題原方程：代入具體數(shù)據(jù)，并對長度歸一化無量綱處34對照標準模型，T就是標準模型中的函數(shù)u，m=0,

a=0,b=1,t0=0,tf=1標準模型：初始條件：零時刻所有位置溫度為30，即u0=30

邊界條件：已知在零位置處任意時間溫度為30，在1位置處，偏導(dǎo)為0pa=u-30,qa=0;pb=0,qb=1

對照標準模型，T就是標準模型中的函數(shù)u，m=0,a=0,b35程序清單functionl51clcclearallglobaluaubalpha=1.0;%cm/sua=30;ub=30;m=0;a=0;b=1;t0=0;tf=1x=linspace(a,b,11);t=linspace(t0,tf,101);sol=pdepe(m,@PDEfun,@ICfun,@BCfun,x,t);u=sol(:,:,1)%surfaceplotofthesolutionfigure;surf(x,t,u);title('Numericalsolutioncomputedwith11meshpoints.');xlabel('Disuancex');ylabel('Timet');ppp=u%solutionprofileatt=1figure;subplot(2,2,1);plot(x,u(1,:),'o-');title('Solutionsatt=0.');xlabel('Distancex');ylabel('u(t=0,)');gridon;subplot(2,2,2);程序清單functionl51u=sol(:,:,136程序清單plot(x,u(50,:),'o-');title('Solutionsatt=0.5.');xlabel('Distancex');ylabel('u(t=0.5)');gridon;subplot(2,2,3);plot(x,u(end,:),'o');title('Solutionsatt=1.');xlabel('Distancex');ylabel('u(t=1,)');gridon;%------------------------------------------------------------------function[c1,f,s]=PDEfun(x,t,u,Du)c1=1f=0.001*Du;s=2*(150-t)-3*Du%------------------------------------------------------------------functionu0=ICfun(x)u0=30;%------------------------------------------------------------------function[pa,qa,pb,qb]=BCfun(xl,ul,xr,ur,t)globaluaubpa=ul-ua;qa=0;pb=0qb=1程序清單plot(x,u(50,:),'o-');c1=37計算結(jié)果計算結(jié)果38分析討論分析討論39Matlab-求解化工常微分方程和偏微分方程課件40Matlab求解化工常微分方程和偏微分方程方利國Matlab求解化工常微分方程和偏微分方程方利國41Matlab求解化工常微分方程和偏微分方程1、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命令1.2初值問題求解1.3邊值問題求解1.4加權(quán)問題求解（自學內(nèi)容）2、偏微分方程（組）求解2.1問題描述及一維動態(tài)PDE方程求解2

.2二維求解Matlab求解化工常微分方程和偏微分方程1、常微分方程（421、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命令

常微分方程：（初值問題）常微分方程：（兩點邊值問題）

1、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命431、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命令

微分方程組：

1、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命441、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命令

高級微分方程：

1、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命451、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命令Matlab調(diào)用命令：ODE45:4-5階龍格庫塔法（非剛性）ODE23:2-3階龍格庫塔法（非剛性）ODE113:可變D-B-M法（非剛性）ODE15S:基于數(shù)值差分的可變階方法法（剛性）ODE23S、ODE23t、ODE23tb（剛性）

1、常微分方程（組）求解1.1問題描述及Matlab調(diào)用命461、常微分方程（組）求解通用調(diào)用格式：[x,y]=ode***(@odefun,xspan,y0,<option,p1,p2..>)X:自變量向量，在實際調(diào)用時取名不一定要用x，也可以用其他名稱，只要前后一致即可。Y:應(yīng)變量向量，在實際調(diào)用時取名不一定要用Y，也可以用其他名稱，只要前后一致即可。***:根據(jù)不同的問題調(diào)用不同格式，如45，23s@odefun:自定義函數(shù)的函數(shù)名，該函數(shù)為Xspan：自變量的積分限，[xa,xb],也可以是離散點，[x0,x1,x2,…xf]y0：應(yīng)變量向量的初值<>:可以沒有該選項，如有，具體應(yīng)用見下面的實際例子

1、常微分方程（組）求解通用調(diào)用格式：471.2初值問題求解例1：已知某高溫物體其溫降過程符合以下規(guī)律，其中溫度T的單位為K，時間

的單位為分鐘，零時刻高溫物體的溫度為2000K，以1分鐘作為時間步長，請計算零時刻以后每隔1分鐘至170分鐘的溫度。單個微分方程1.2初值問題求解例1：已知某高溫物體其溫降過程符合以下48functionxODEs%鐵球從2000K降溫曲線,在7.0版本上調(diào)試通過%由華南理工大學方利國編寫，2012年2月29日%歡迎讀者調(diào)用，如有問題請告知lgfang@clearall；clcy0=[2000];[x1,y1]=ode45(@f,[0:1:170],y0);%0到170分鐘，每分鐘一個計算點[x2,y2]=ode23(@f,[0:1:170],y0);plot(x1,y1,'r-')xlabel('時間，M')ylabel('溫度，K')holdondisp('Resultsbyusingode45()：')disp('xy(1)')disp([x1y1])disp('Resultsbyusingode23()：')disp('xy(2)')disp([x2y2])plot(x2,y2,'k:')%------------------------------------------------------------------functiondy=f(x,y)%定義降溫速率的微分方程%dy=0.04*y(1)-100;dy=-0.04*exp(0.001*(y(1)-300))*(-300+y(1));functionxODEs49Resultsbyusingode45()：

xy(1)1.0e+003*Columns1through1300.01000.02000.03000.04000.05000.06000.07000.08000.09000.10000.11000.12002.00000.87880.61330.48800.41810.37620.34980.33280.32180.31450.30970.30650.3043Columns14through180.13000.14000.15000.16000.17000.30290.30190.30130.30090.3006Resultsbyusingode23()：

xy(2)1.0e+003*Columns1through1300.01000.02000.03000.04000.05000.06000.07000.08000.09000.10000.11000.12002.00000.87790.61240.48730.41760.37560.34930.33230.32130.31400.30920.30610.3041Columns14through180.13000.14000.15000.16000.17000.30270.30180.30120.30080.3005Resultsbyusingode45()：50Matlab-求解化工常微分方程和偏微分方程課件511.2初值問題求解該問題相當與一個自變量，兩個應(yīng)變量問題，已知初值及微分表達式，可以利用ODE45求解。微分方程組求解1.2初值問題求解該問題相當與一個自變量，兩個應(yīng)變量問題52程序代碼functionuvDEs%微生物消亡問題計算,在7.0版本上調(diào)試通過%由華南理工大學方利國編寫，2012年3月12日%歡迎讀者調(diào)用，如有問題請告知lgfang@clearall;Clcy0=[1.61.2];[x1,y1]=ode45(@f,[0:0.1:10],y0);%0到3分鐘，每0.1分鐘一個計算點u=y1(:,1);v=y1(:,2);plot(x1,u,'r-')xlabel('時間，M')ylabel('微生物濃度')holdonplot(x1,v,'k:')disp('Resultsbyusingode45()：')disp('xuv')disp([x1y1])%------------------------------------------------------------------functiondy=f(x,y)%定義降溫速率的微分方程f1=0.09*y(1)*(1-y(1)/20)-0.45*y(1)*y(2);f2=0.06*y(2)*(1-y(2)/15)-0.001*y(1)*y(2);dy=[f1;f2];程序代碼functionuvDEs531.2初值問題求解

例3：當X較大時，兩種方法計算結(jié)果有較大不同，為什么？單個微分方程有零點問題？1.2初值問題求解例3：當X較大時，兩種方法計算結(jié)果有54functionL43ODEs%在7.0版本上調(diào)試通過%由華南理工大學方利國編寫，2012年2月29日%歡迎讀者調(diào)用，如有問題請告知lgfang@clearallclcy0=[1];[x1,y1]=ode45(@f,[0:0.05:10],y0);%0到10，每0.05間隔一個計算點[x2,y2]=ode23(@f,[0:0.05:10],y0);%0到10，每0.05間隔一個計算點plot(x1,y1,'r-')xlabel('x')ylabel('y')holdondisp('Resultsbyusingode45()：')disp('xy(1)')disp([x1y1])disp('Resultsbyusingode23()：')disp('xy(2)')disp([x2y2])plot(x2,y2,'b--')%------------------------------------------------------------------functiondy=f(x,y)%定義微分方程dy=y^2*cos(x);functionL43ODEs55計算值xy(1)Columns1through1300.05000.10000.15000.20000.25000.30000.35000.40000.45000.50000.55000.60001.00001.05261.11091.17571.24791.32871.41951.52181.63781.76981.92102.09512.2970Columns14through260.65000.70000.75000.80000.85000.90000.95001.00001.05001.10001.15001.20001.25002.53292.81073.14113.53814.02064.61525.35986.30787.54359.192811.461414.728319.5991Columns27through291.30001.35001.400027.428341.203068.6630計算值xy(1)56高階微分方程求解求解思路：將高階微分方程通過變量轉(zhuǎn)換，轉(zhuǎn)變成一級微分方程組進行求解。例4：高階微分方程求解求解思路：將高階微分方程通過變量轉(zhuǎn)換，轉(zhuǎn)變成57高階微分方程求解程序及解Columns1through1300.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.10001.200000.10020.20150.30520.41290.52630.64760.77890.92301.08271.26141.46271.69081.00001.00531.02271.05441.10261.16981.25881.37231.51361.68611.89352.13982.4296Columns14through211.30001.40001.50001.60001.70001.80001.90002.00001.95022.24612.58412.97053.41233.91714.49365.15132.76773.15953.61104.12864.71965.39186.15417.0159高階微分方程求解程序及解Columns1through58剛性方程求解有些微分組的系數(shù)變化很大，這時用ODE45就很難收斂求解，這時可用專門解決此類微分方程的ODE23S來求解，需要注意的是在解的圖像繪制時，也需要考慮數(shù)值的波動幅度很大，需要引入對數(shù)坐標。例5：

dy1/dx=-0.03*y1+1e4*y2*y3dy2/dx=0.03*y1-2e4*y2*y3-5e7*y2^2dy3/dx=5e7*y2^2y1(0)=1,y2(0)=0,y3(0)=0剛性方程求解有些微分組的系數(shù)變化很大，這時用ODE45就很難59functiongangxinDEs%剛性問題計算,在7.0版本上調(diào)試通過%由華南理工大學方利國編寫，2012年3月13日%歡迎讀者調(diào)用，如有問題請告知lgfang@%dy1=-0.03*y1+1e4*y2*y3%dy2=0.03*y1-2e4*y2*y3-5e7*y2^2%dy3=5e7*y2^2clearallclcfigurexspan=[06*logspace(-6,6)];y0=[100];[x1,y1]=ode15s(@f,xspan,y0);%用ode45計算剛性方程，可能有問題u=y1(:,1);v=1e4*y1(:,2);w=y1(:,3);semilogx(x1,u,'r-','linewidth',2)xlabel('x')ylabel('1e4*v')holdonsemilogx(x1,v,'k:','linewidth',2)holdonsemilogx(x1,w,'g-','linewidth',2)gridaxis([10^(-10)10^10-0.21.2])legend('u','v','w')disp('Resultsbyusingode45()：')disp('xuvw')formatlongdisp([x1y1])%------------------------------------------------------------------functiondy=f(x,y)%定義微分方程f1=-0.03*y(1)+1e4*y(2)*y(3);f2=0.03*y(1)-2e4*y(2)*y(3)-5e7*y(2)^2;f3=5e7*y(2)^2;dy=[f1;f2;f3];functiongangxinDEsdisp('x601.3邊值問題求解邊值問題相對于初值問題而言，多了一個端點的約束，如果在高階或微分方程組中端點約束過多，微分方程組可能無解，端點約束有一定限制。可以通過建立離散的方程組，再利用ODE45進行求解，但可以利用MATLAB的專用工具求解最好。下面介紹ODE-BVPs的求解器，主要有bvpinit,bvp4c,deval,solinit等。通過實際例子介紹這些內(nèi)部函數(shù)的功能。1.3邊值問題求解邊值問題相對于初值問題而言，多了一個端點611.3邊值問題求解solinit=bvpinit（x,yinit):產(chǎn)生在初始網(wǎng)格上的初始解，以便bvp4c調(diào)用，其中x為自變量網(wǎng)格，yinit為對應(yīng)函數(shù)的初值。sol=bvp4c(@odefun,@BCfun,solinit,<option,p1,p2..>)@Bcfun:為定義邊界條件方程，Bcfun(ya,yb),其中ya、yb分別表示左右邊界。其他符號意義同上。deval(sol,xint):計算任意點處的函數(shù)值。1.3邊值問題求解solinit=bvpinit（x,yi621.3邊值問題求解

例6：1.3邊值問題求解例6：63程序functionBVP4c1%求解兩點邊值問題的示例在7.0版本上調(diào)試通過%由華南理工大學方利國編寫，2012年3月13日%歡迎讀者調(diào)用，如有問題請告知lgfang@clearallclca=0;b=10;solinit=bvpinit(linspace(a,b,101),[00]);sol=bvp4c(@ODEfun,@BCfun,solinit);formatshortx=[0:0.1:10];y=deval(sol,x);y1=y(1,:)y2=y(2,:)plot(x,y1,'r-')xlabel('x')ylabel('y')holdongridplot(x,y2,'k:')legend('y','dy')disp('Resultsbyusingbvp4c：')disp('xydy')disp([x1y1])%------------------------------------------------------------------functiondy=ODEfun(x,y)f1=y(2);f2=0.05*(1+x^2)*y(1)+2;dy=[f1;f2];%------------------------------------------------------------------functionbc=BCfun(ya,yb)bc=[ya(1)-40;yb(1)-80];Resultsbyusingbvp4c：

xydyColumns1through1300.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.10001.200039.999838.172236.384034.634732.924431.253129.621428.029826.479124.970223.503822.081020.7025-18.4739-18.0779-17.6872-17.2985-16.9088-16.5158-16.1176-15.7125-15.2996-14.8781-14.4476-14.0080-13.5597程序functionBVP4c1disp('Results64計算結(jié)果計算結(jié)果651.3邊值問題求解

例7：1.3邊值問題求解例7：66主要程序

functiondy=ODEfun(x,y)f1=y(2);f2=4*y(1);dy=[f1;f2];%------------------------------------------------------------------functionbc=BCfun(ya,yb)bc=[ya(1)-1;yb(1)-exp(4)];主要程序

functiondy=ODEfun(x,y)67functionBVP4c2%求解兩點邊值問題的示例在7.0版本上調(diào)試通過%由華南理工大學方利國編寫，2012年3月13日%歡迎讀者調(diào)用，如有問題請告知lgfang@clearallclca=0;b=2solinit=bvpinit(linspace(a,b,21),[00]);sol=bvp4c(@ODEfun,@BCfun,solinit);formatshortx=[0:0.1:2];y=deval(sol,x);y1=y(1,:);y2=y(2,:);y3=exp(2*x);plot(x,y1,'r-','linewidth',2)xlabel('x')ylabel('y')holdongridplot(x,y2,'k:','linewidth',2)holdonplot(x,y3,'b:','linewidth',2)legend('y','dy','real-y')disp('Resultsbyusingbvp4c：')disp('xydy')disp([x;y1;y2])%------------------------------------------------------------------functiondy=ODEfun(x,y)f1=y(2);f2=4*y(1);dy=[f1;f2];%-----