Matlab求解边值问题方法例题.pdf

Matlab求解边值问题方法 ：bvp4c 函数 求解边值问题方法 p 函数 1. 把待解的问题转化为标准边值问题 (,) () ,() )0 yfx y gyayb 2. 因为边值问题可以多解，所以需要为期望解指定一个初始 猜测解。该猜测解网（Mesh）包括区间a, b内的一组网 点（Mesh points）和网点上的解S(x) 3 根据原微分方程构造残差函数 () () ( () ) SfS 3. 根据原微分方程构造残差函数 4. 利用原微分方程和边界条件，借助迭代不断产生新的S(x) ， () () (, () ) rx S x f x S x 使残差不断减小，从而获得满足精度要求的解Matlab求解边值问题的基本指令 求解边值问题的基本指令 solinit=bvpinit(x,v,parameters) 生成b4 调用指令所必须的 “ 解猜测网 ” 生成bvp4c调用指令所必须的 “ 解猜测网 ” sol=bvp4c(odefun bcfun solinit options p1 p2 ) sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,) 给出微分方程边值问题的近似解 sxint=deval(sol ，xint) 计算微分方程积分区间内任何一点的解值初始解生成函数 ：bvpinit() 初始解生成函数 p( ) solinit=bvpinit(x,v,parameters) x指定边界区间a,b上的初始网络，通常是等距排列的（1 M）一维数组。 注意：使x(1)=a ，x(end)=b；格点要单调排列。 是对解的初始猜测 v是对解的初始猜测 solinit（可以取别的任意名）是“解猜测网（Mesh）”。 它是一个结构体，带如下两个域： 是表示初始网格有序节点的 维数组 并 solinit.x是表示初始网格有序节点的 （1 M）一维数组 ， 并 且 solinit.x(1)一定是a ，solinit.x(end)一定是b 。M不宜取得太大，10 数量级左右即可。 solinit.y是表示网点上微分方程解的猜测值的（N M）二维数组。 solinit.y(:,i)表示节点solinit.x(i)处的解的猜测值。边值问题求解指令 ：bvp4c() sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,) 边值问题求解指令 p () p ( , , , p ,p ,p , ) 输入参数： odefun是计算导数的m函数文件。该函数的基本形式为： dydx=odefun(x y parameters p1 p2 ) 在此 自变量x是标量 y dydx=odefun(x,y,parameters,p1,p2,) ，在此 ， 自变量x是标量 ，y ， dydx是列向量。 bcfun是计算边界条件下残数的m函数文件。其基本形式为： bf ( b t 12 ) 文件输入宗量 b是边界 res=bcfun(ya,yb,parameters,p1,p2,) ， 文件输入宗量ya,yb是边界 条件列向量，分别代表y 在a 和b处的值。res是边界条件满足时的残数 列向量。注意：例如odefun函数的输入宗量中包含若干“未知”和 “已知”参数，那么不管在边界条件计算中是否用到，它们都应作为 bcfun的输入宗量。 输入宗量options是用来改变bvp4c算法的控制参数的。在最基本用法 p p 中，它可以缺省，此时一般可以获得比较满意的边值问题解。如需更 改可采用bvpset函数，使用方法同odeset函数。 输入宗量p1,p2等表示希望向被解微分方程传递的已知参数 。 如果无须 输入宗量p1,p2等表示希望向被解微分方程传递的已知参数 。 如果无须 向微分方程传递参数，它们可以缺省。边值问题求解指令 ：bvp4c() 边值问题求解指令 p () 输出参数： 输出变量sol是一个结构体 sol.x是指令bvp4c所采用的网格节点； sol.y 是y(x) 在sol.x网点上的近似解值； sol yp 是y(x) 在sol x网点上的近似解值 sol.yp 是y(x) 在sol.x网点上的近似解值 ； sol.parameters是微分方程所包含的未知参数的近似解值。 当被解微分方程包含未知参数时 ， 该域存在 。 当被解微分方程包含未知参数时 ， 该域存在 。边值问题的求解 边值问题的求解 10 (0) (1) 0 yy yy 求解两点边值问题： 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 令： 原方程组等价于以下标准形式的 方程组： 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函数 function dydx=ODEfun(x y) 12 21 10 yy yy function dydx=ODEfun(x,y) dydx=y(2);y(1)+10; % 定义BCfun函数 function bc=BCfun(ya,yb) (0) 0 (1) 0 yy 边界条件为： bc=ya(1);yb(1); 11 (0) 0, (1) 0 yy边值问题的求解 2 边值问题的求解 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 令： 原方程组等价于以下标准形式的 方程组： 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函数 function dydx=ODEfun(x,y) 12 2 21 (1 ) 1 yy yx y dydx=y(2);(1+x2)*y(1)+1; % 定义BCfun函数 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边值问题的求解 边值问题的求解 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( 解曲线,初始网格点解) % 定义ODEfun函数 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函数 function bc=BCfun(ya,yb,c) 121 令 bc=ya(1);yb(1)+2;边值问题的求解 solinit=bvpinit(linspace(0,pi,10),1;1,lmb); 边值问题的求解 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(解曲线,初始网格点解) 边界条件： q 0 ) ( 0 ) 0 ( 1 ) 0 ( z z z g( 解曲线 , 初始网格点解 ) % 定义ODEfun函数 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函数 function bc=BCfun(ya,yb,lmb) bc=ya(1)-1;ya(2);yb(2); y ( ) ;y ( );y ( ); 本例中，微分方程与参数 的数值有关。一般而言，对于任意的 值，该问题无解， 但对于特殊的 值（特征值），它存在一个解，这也称为微分方程的特征值问题。 对于此问题，可在bvpinit中提供参数的猜测值，然后重复求解BVP 得到所需的 参数，返回参数为sol.parameters边值问题的求解 infinity=6; solinit=bvpinit(linspace(0 infinity 5) 0 0 1); 如果取1，计算结果如何? 边值问题的求解 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/deta),shg % 定义ODEfun函数 边界条件： (0) 0, (0) 0, () 1 , ff f % 定义ODEfun函数 function dfdeta=ODEfun(eta,f) beta=0.5; dfdeta=f(2);f(3);-f(1)*f(3)-beta*(1-f(2)2); % 定义BCfun函数 function res=BCfun(f0,finf) res=f0(1);f0(2);finf(2)-1;工程应用中我们经常遇到边值问题，这些是那些ode*函数无能为力的， 当然我们可以自己编写函数求解(shooting) ，但是不是能力所及的. 中提供了 解算器 Matlab中提供了 bvp解算器 。 solinit = bvpinit(x, yinit, params) sol = bvpsolver(odefun,bcfun,solinit,options) sol bvpsolver(odefun,bcfun,solinit,options) 由于边值问题可能有多解，为了便于我们确定那个解是我们需要的，所 以必须使用bvpinit函数对初值进行估计 解算器(bl ) 解算器(bvpsolver) ： Matlab中提供了bvp4c 和bvp5c，后者误差控制更好些 输 入 参数 ： 输参 数 x：需要计算的网格点，相当于ode* 的tspan yinit：猜测的值，可以是具体值，也可以是函数，类似与 ode* 的 x0 params 其它未知参数 也是 个猜测值 params ： 其它未知参数 ， 也是 一 个猜测值 odefun：描述边值问题微分方程的函数句柄 bcfun：边值函数，一般是双边值(x 的上下限即认为两个边界)，支持多 ( ) 边值 solinit ：由bvpinit 生成的初始化网格 options ：BVP解算器优化参数 可以通过bvpset设置 options ：BVP解算器优化参数 ， 可以通过bvpset设置 输出参数：大部分同理ode45