Matlab求解边值问题方法+例题.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,把待解的问题转化为标准边值问题,因为边值问题可以多解，所以需要为期望解指定一个初始猜测解。该猜测解网（,Mesh,）包括区间,a,b,内的一组网点（,Mesh points,）和网点上的解,S(x),根据原微分方程构造残差函数,利用原微分方程和边界条件，借助迭代不断产生新的,S(x),，使残差不断减小，从而获得满足精度要求的解,Matlab,求解边值问题方法：,bvp4c,函数,solinit=bvpinit(x,v,parameters),生成,bvp4c,调用指令所必须的“解猜测网”,sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,),给出微分方程边值问题的近似解,sxint=deval(sol,，,xint),计算微分方程积分区间内任何一点的解值,Matlab求解边值问题,的基本指令,solinit=bvpinit(x,v,parameters),x,指定边界区间,a,b,上的初始网络，通常是等距排列的（,1M,）一维数组。,注意：使,x(1)=a,，,x(end)=b,；格点要单调排列。,v,是对解的初始猜测,solinit,（可以取别的任意名）是“解猜测网（,Mesh,）”。,它是一个结构体，带如下两个域：,solinit.x,是表示初始网格有序节点的（,1M,）一维数组，并且,solinit.x(1),一定是,a,，,solinit.x(end),一定是,b,。,M,不宜取得太大，,10,数量级左右即可。,solinit.y,是表示网点上微分方程解的猜测值的（,NM,）二维数组。,solinit.y(:,i),表示节点,solinit.x(i),处的解的猜测值。,初始解生成函数：,bvpinit(),输出参数：,输出变量,sol,是一个结构体,sol.x,是指令,bvp4c,所采用的网格节点；,sol.y,是,y(x),在,sol.x,网点上的近似解值；,sol.yp,是,y(x),在,sol.x,网点上的近似解值；,sol.parameters,是微分方程所包含的未知参数的近似解值。,当被解微分方程包含未知参数时，该域存在。,边值问题求解指令：,bvp4c,(),原方程组等价于以下标准形式的方程组：,solinit=bvpinit(linspace(0,1,10),1 0);,sol=bvp4c(ODEfun,BCfun,solinit);,x=0:0.05:0.5;,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),dydx=y(2);y(1)+10;,%,定义,BCfun,函数,function bc=BCfun(ya,yb),bc=ya(1);yb(1);,求解两点边值问题：,令：,边界条件为：,边值问题的求解,原方程组等价于以下标准形式的方程组：,solinit=bvpinit(linspace(0,1,10),0 1);,sol=bvp4c(ODEfun,BCfun,solinit);,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.,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),legend boxoff,%,定义,ODEfun,函数,function dydx=ODEfun(x,y),dydx=y(2);(1+x2)*y(1)+1;,%,定义,BCfun,函数,function bc=BCfun(ya,yb),bc=ya(1)-1;yb(1)-3;,求解：,令：,边界条件为：,边值问题的求解,c=1;,solinit=bvpinit(linspace(0,4,10),1 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),legend(,解曲线,初始网格点解,),%,定义,ODEfun,函数,function dydx=ODEfun(x,y,c),dydx=y(2);-c*abs(y(1);,%,定义,BCfun,函数,function bc=BCfun(ya,yb,c),bc=ya(1);yb(1)+2;,求解：,令：,边值问题的求解,infinity=6;,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&Keller report f(0)=0.92768.n),fprintf(Value computed here is f(0)=%7.5f.n,f(3,1),clf reset,plot(eta,f(2,:);,axis(0 infinity 0 1.4);,title.,(Falkner-Skan equation,positive wall shear,beta=0.5.),xlabel(eta),ylabel(df/deta),shg,%,定义,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;,求解：,边界条件：,如果取,1,，计算结果如何,?,边值问题的求解,