Matlab求解邊值問(wèn)題方法例題.pdf

《Matlab求解邊值問(wèn)題方法例題.pdf》由會(huì)員分享，可在線(xiàn)閱讀，更多相關(guān)《Matlab求解邊值問(wèn)題方法例題.pdf（11頁(yè)珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

Matlab求解邊值問(wèn)題方法 ：bvp4c 函數(shù) 求解邊值問(wèn)題方法 p 函數(shù) 1. 把待解的問(wèn)題轉(zhuǎn)化為標(biāo)準(zhǔn)邊值問(wèn)題 (,) (() ,() )0 ? ? ? ? ? ? yfx y gyayb 2. 因?yàn)檫呏祮?wèn)題可以多解，所以需要為期望解指定一個(gè)初始 猜測(cè)解。該猜測(cè)解網(wǎng)（Mesh）包括區(qū)間[a, b]內(nèi)的一組網(wǎng) [] 點(diǎn)（Mesh points）和網(wǎng)點(diǎn)上的解S(x) 3 根據(jù)原微分方程構(gòu)造殘差函數(shù) () () ( () ) ? SfS 3. 根據(jù)原微分方程構(gòu)造殘差函數(shù) 4. 利用原微分方程和邊界條件，借助迭代不斷產(chǎn)生新的S(x) ， () () (, () ) ? ?? rx S x f x S x 使殘差不斷減小，從而獲得滿(mǎn)足精度要求的解Matlab求解邊值問(wèn)題的基本指令 求解邊值問(wèn)題的基本指令 solinit=bvpinit(x,v,parameters) 生成b4 調(diào)用指令所必須的 “ 解猜測(cè)網(wǎng) ” 生成bvp4c調(diào)用指令所必須的 “ 解猜測(cè)網(wǎng) ” sol=bvp4c(odefun bcfun solinit options p1 p2 ) sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) 給出微分方程邊值問(wèn)題的近似解 sxint=deval(sol ，xint) 計(jì)算微分方程積分區(qū)間內(nèi)任何一點(diǎn)的解值初始解生成函數(shù) ：bvpinit() 初始解生成函數(shù) p( ) solinit=bvpinit(x,v,parameters) ? x指定邊界區(qū)間[a,b]上的初始網(wǎng)絡(luò)，通常是等距排列的（1 M）一維數(shù)組。 注意：使x(1)=a ，x(end)=b；格點(diǎn)要單調(diào)排列。 ? 是對(duì)解的初始猜測(cè) ? v是對(duì)解的初始猜測(cè) ? solinit（可以取別的任意名）是“解猜測(cè)網(wǎng)（Mesh）”。 它是一個(gè)結(jié)構(gòu)體，帶如下兩個(gè)域： 是表示初始網(wǎng)格有序節(jié)點(diǎn)的 維數(shù)組 并 ? solinit.x是表示初始網(wǎng)格有序節(jié)點(diǎn)的 （1 M）一維數(shù)組 ， 并 且 solinit.x(1)一定是a ，solinit.x(end)一定是b 。M不宜取得太大，10 數(shù)量級(jí)左右即可。 ? solinit.y是表示網(wǎng)點(diǎn)上微分方程解的猜測(cè)值的（N M）二維數(shù)組。 solinit.y(:,i)表示節(jié)點(diǎn)solinit.x(i)處的解的猜測(cè)值。邊值問(wèn)題求解指令 ：bvp4c() sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) 邊值問(wèn)題求解指令 p () p ( , , , p ,p ,p , ) 輸入?yún)?shù)： ? odefun是計(jì)算導(dǎo)數(shù)的m函數(shù)文件。該函數(shù)的基本形式為： dydx=odefun(x y parameters p1 p2 ) 在此 自變量x是標(biāo)量 y dydx=odefun(x,y,parameters,p1,p2,…) ，在此 ， 自變量x是標(biāo)量 ，y ， dydx是列向量。 ? bcfun是計(jì)算邊界條件下殘數(shù)的m函數(shù)文件。其基本形式為： bf ( b t 12 ) 文件輸入宗量 b是邊界 res=bcfun(ya,yb,parameters,p1,p2,…) ， 文件輸入宗量ya,yb是邊界 條件列向量，分別代表y 在a 和b處的值。res是邊界條件滿(mǎn)足時(shí)的殘數(shù) 列向量。注意：例如odefun函數(shù)的輸入宗量中包含若干“未知”和 “已知”參數(shù)，那么不管在邊界條件計(jì)算中是否用到，它們都應(yīng)作為 bcfun的輸入宗量。 ? 輸入宗量options是用來(lái)改變bvp4c算法的控制參數(shù)的。在最基本用法 p p 中，它可以缺省，此時(shí)一般可以獲得比較滿(mǎn)意的邊值問(wèn)題解。如需更 改可采用bvpset函數(shù)，使用方法同odeset函數(shù)。 ? 輸入宗量p1,p2等表示希望向被解微分方程傳遞的已知參數(shù) 。 如果無(wú)須 ? 輸入宗量p1,p2等表示希望向被解微分方程傳遞的已知參數(shù) 。 如果無(wú)須 向微分方程傳遞參數(shù)，它們可以缺省。邊值問(wèn)題求解指令 ：bvp4c() 邊值問(wèn)題求解指令 p () 輸出參數(shù)： 輸出變量sol是一個(gè)結(jié)構(gòu)體 ?sol.x是指令bvp4c所采用的網(wǎng)格節(jié)點(diǎn)； ?sol.y 是y(x) 在sol.x網(wǎng)點(diǎn)上的近似解值； ?sol yp 是y(x) 在sol x網(wǎng)點(diǎn)上的近似解值 ?sol.yp 是y(x) 在sol.x網(wǎng)點(diǎn)上的近似解值 ； ?sol.parameters是微分方程所包含的未知參數(shù)的近似解值。 當(dāng)被解微分方程包含未知參數(shù)時(shí) ， 該域存在 。 當(dāng)被解微分方程包含未知參數(shù)時(shí) ， 該域存在 。邊值問(wèn)題的求解 ?? ? 邊值問(wèn)題的求解 10 (0) (1) 0 ?? ? ? ? ? ? ? ? yy yy 求解兩點(diǎn)邊值問(wèn)題： solinit=bvpinit(linspace(0,1,10),[1 0]); sol=bvp4c(@ODEfun,@BCfun,solinit); x=[0:0.05:0.5]; y=deval(sol x); 121 , ? ?? yyyy 令： 原方程組等價(jià)于以下標(biāo)準(zhǔn)形式的 方程組： y=deval(sol,x); xP=[0:0.1:0.5]; yP=[0 -0.41286057 -0.729740656... -0.95385538 -1.08743325 -1.13181116]; ? ? ? ] plot(xP,yP,o,x,y(1,:),r-) legend(Analytical Solution,Numerical Solution) % 定義ODEfun函數(shù) function dydx=ODEfun(x y) 12 21 10 ? ? ? ? ? ? ?? ? ? yy yy function dydx=ODEfun(x,y) dydx=[y(2);y(1)+10]; % 定義BCfun函數(shù) function bc=BCfun(ya,yb) (0) 0 (1) 0 ?? yy 邊界條件為： bc=[ya(1);yb(1)]; 11 (0) 0, (1) 0 ?? yy邊值問(wèn)題的求解 2 ? 邊值問(wèn)題的求解 solinit=bvpinit(linspace(0,1,10),[0 1]); sol=bvp4c(@ODEfun,@BCfun,solinit); x=[0:0 1:1]; 求解： 2 (1 ) 1 (0) 1, (1) 3 ?? ? ?? ? ? ?? ? y x y yy x=[0:0.1:1]; y=deval(sol,x); xP=[0:0.1:1.0]; yP=[1 1.0743 1.1695 1.2869 1.4284... 121 , ? ?? yyyy 令： 原方程組等價(jià)于以下標(biāo)準(zhǔn)形式的 方程組： y[ 1.5965 1.7947 2.0274 2.3004 2.6214 3]; plot(xP,yP,o,x,y(1,:),r-) legend(Analytical Solution,Numerical Solution,... location Northwest) ? ? location , Northwest ) legend boxoff % 定義ODEfun函數(shù) function dydx=ODEfun(x,y) 12 2 21 (1 ) 1 ? ? ? ? ? ? ?? ? ? ? yy yx y dydx=[y(2);(1+x^2)*y(1)+1]; % 定義BCfun函數(shù) function bc=BCfun(ya,yb) bc=[ya(1) 1;yb(1) 3]; 邊界條件為： ? ? ? ? 0 1 ) 0 ( 1 y bc=[ya(1)-1;yb(1)-3]; ? ? ? ? 0 3 ) 1 ( ) ( 1 1 y y邊值問(wèn)題的求解 邊值問(wèn)題的求解 c=1; solinit=bvpinit(linspace(0,4,10),[1 1]); sol=bvp4c(@ODEfun @BCfun solinit [] c); 求解： 0 ?? ? 1 sol bvp4c(@ODEfun,@BCfun,solinit,[],c); x=[0:0.1:4]; y=deval(sol,x); plot(x,y(1,:),b-,sol.x,sol.y(1,:),ro) 0 ? ?? z cz 2 ) 4 ( , 0 ) 0 ( ? ? ? z z 1 ? c plot(x,y(1,:), b ,sol.x,sol.y(1,:), ro ) legend( 解曲線(xiàn),初始網(wǎng)格點(diǎn)解) % 定義ODEfun函數(shù) function dydx=ODEfun(x,y,c) 121 , ? ?? z zz z 令 ： 2 ) 4 ( , 0 ) 0 ( z z y( , y , ) dydx=[y(2);-c*abs(y(1))]; % 定義BCfun函數(shù) function bc=BCfun(ya,yb,c) 121 令 bc=[ya(1);yb(1)+2];邊值問(wèn)題的求解 solinit=bvpinit(linspace(0,pi,10),[1;1],lmb); 邊值問(wèn)題的求解 p(p( , p ,) , [ ; ] ,) ; opts=bvpset(Stats,on); sol=bvp4c(@ODEfun,@BCfun,solinit,opts); lambda=sol.parameters [0 i/60 i] 求解： (2c o s 2 )0 ?? ?? ? zqx z ? 15 ? q x=[0:pi/60:pi]; y=deval(sol,x); plot(x,y(1,:),b-,sol.x,sol.y(1,:),ro) legend(解曲線(xiàn),初始網(wǎng)格點(diǎn)解) 邊界條件： q 0 ) ( 0 ) 0 ( 1 ) 0 ( ? ? ? z z z g( 解曲線(xiàn) , 初始網(wǎng)格點(diǎn)解 ) % 定義ODEfun函數(shù) function dydx=ODEfun(x,y,lmb) q=15; dydx=[y(2); (lmb 2*q*cos(2*x))*y(1)]; 0 ) ( , 0 ) 0 ( , 1 ) 0 ( ? ? ? ? z z z dydx=[y(2);-(lmb-2*q*cos(2*x))*y(1)]; % 定義BCfun函數(shù) function bc=BCfun(ya,yb,lmb) bc=[ya(1)-1;ya(2);yb(2)]; [y ( ) ;y ( );y ( )]; 本例中，微分方程與參數(shù) λ 的數(shù)值有關(guān)。一般而言，對(duì)于任意的 λ 值，該問(wèn)題無(wú)解， 但對(duì)于特殊的 λ值（特征值），它存在一個(gè)解，這也稱(chēng)為微分方程的特征值問(wèn)題。 對(duì)于此問(wèn)題，可在bvpinit中提供參數(shù)的猜測(cè)值，然后重復(fù)求解BVP 得到所需的 參數(shù)，返回參數(shù)為sol.parameters邊值問(wèn)題的求解 infinity=6; solinit=bvpinit(linspace(0 infinity 5) [0 0 1]); 如果取1，計(jì)算結(jié)果如何? 邊值問(wèn)題的求解 solinit=bvpinit(linspace(0,infinity,5),[0 0 1]); options=bvpset(stats,on); sol=bvp4c(@ODEfun,@BCfun,solinit,options); eta=sol.x; f=sol.y; fprintf(\n); fprintf(Cebeci axis([0 infinity 0 1.4]); 0.5 ? ? title... (Falkner-Skan equation, positive wall shear,\beta=0.5.) xlabel(\eta),ylabel(df/d\eta),shg % 定義ODEfun函數(shù) 邊界條件： (0) 0, (0) 0, () 1 , ? ?? ? ?? ? ff f?? % 定義ODEfun函數(shù) function dfdeta=ODEfun(eta,f) beta=0.5; dfdeta=[f(2);f(3);-f(1)*f(3)-beta*(1-f(2)^2)]; % 定義BCfun函數(shù) function res=BCfun(f0,finf) res=[f0(1);f0(2);finf(2)-1];工程應(yīng)用中我們經(jīng)常遇到邊值問(wèn)題，這些是那些ode**函數(shù)無(wú)能為力的， 當(dāng)然我們可以自己編寫(xiě)函數(shù)求解(shooting) ，但是不是能力所及的. 中提供了 解算器 Matlab中提供了 bvp解算器 。 solinit = bvpinit(x, yinit, params) sol = bvpsolver(odefun,bcfun,solinit,options) sol bvpsolver(odefun,bcfun,solinit,options) 由于邊值問(wèn)題可能有多解，為了便于我們確定那個(gè)解是我們需要的，所 以必須使用bvpinit函數(shù)對(duì)初值進(jìn)行估計(jì) 解算器(bl ) 解算器(bvpsolver) ： Matlab中提供了bvp4c 和bvp5c，后者誤差控制更好些 輸 入 參數(shù) ： 輸參 數(shù) x：需要計(jì)算的網(wǎng)格點(diǎn)，相當(dāng)于ode** 的tspan yinit：猜測(cè)的值，可以是具體值，也可以是函數(shù)，類(lèi)似與 ode** 的 x0 params 其它未知參數(shù) 也是 個(gè)猜測(cè)值 params ： 其它未知參數(shù) ， 也是 一 個(gè)猜測(cè)值 odefun：描述邊值問(wèn)題微分方程的函數(shù)句柄 bcfun：邊值函數(shù)，一般是雙邊值(x 的上下限即認(rèn)為兩個(gè)邊界)，支持多 ( ) 邊值 solinit ：由bvpinit 生成的初始化網(wǎng)格 options ：BVP解算器優(yōu)化參數(shù) 可以通過(guò)bvpset設(shè)置 options ：BVP解算器優(yōu)化參數(shù) ， 可以通過(guò)bvpset設(shè)置 輸出參數(shù)：大部分同理ode45