工程数学作业第十二次方健(1213)

姓名：方健学号：652081701073专业：化学工程班级：化工1503问题：利用MATLAB对于边值问题求解器bvp4c求解一个常微分方程组，并且对于bvp4c函数的用法经行熟悉与应用。解决两点边值问题(BVPs)常微分方程例如求解如下的问题：''22101,13yxyyy令'121,yyyy所以原方程等价于以下标准形式的方程组：'12'22121yyyxy边界条件：11(0)10(1)30yy根据Matlab求解边值问题的基本指令solinit＝bvpinit（x,v,parameters）生成bvp4c调用指令所必须的“解猜测网，利用matlab中的help知道sol＝bvp4c（odefun，bcfun，solinit，options，p1，p2，…）给出微分方程边值问题的近似解sxint＝deval（sol，xint）计算微分方程积分区间内任何一点的解值。sol＝bvp4c（odefun，bcfun，solinit，options，p1，p2，…）输入参数：odefun是计算导数的M函数文件。该函数的基本形式为：dydx＝odefun(x,y,parameters,p1,p2，…)，在此，自变量x是标量，y，dydx是列向量。其基本形式与初值问题的odefile形式相同。bcfun是计算边界条件下残数的M函数文件。其基本形式为：res＝bcfun(ya,yb,parameters,p1,p2,……)，文件输入宗量ya,yb是边界条件列向量,分别代表y在a和b处的值。res是边界条件满足时的残数列向量。这里有一个需要注意：例如odefun函数的输入宗量中包含若干“未知”和“已知”参数，那么不管在边界条件计算中是否用到，它们都应作为bcfun的输入宗量。根据以上分析可以编写以下命令：函数BCfun.m文件：functionbc=BCfun(ya,yb)bc=[ya(1)-1;yb(1)-3];ODEfun.m文件：functiondydx=ODEfun(x,y)dydx=[y(2);(2+x^2)*y(1)+1];程序文件b_bvp.m文件:functionb_bvpclearallclcsolinit=bvpinit(linspace(0,1,10),[01]);sol=bvp4c(@ODEfun,@BCfun,solinit);x=[0:0.1:1];y=deval(sol,x);yP=[11.07431.16951.28691.4284...1.59651.79472.02742.30042.62143];xP=[0:0.1:1.0];plot(xP,yP,'o',x,y(1,:),'r-')legend('AnalyticalSolution','NumericalSolution','location','Northwest')legendboxoff于是在commandwindow中运行b_bvpy=Columns1through102.00001.86921.78601.74941.75941.81761.92682.09172.31882.6175-1.5537-1.0670-0.5980-0.13410.33760.83081.36111.94702.61123.3820Column113.00004.2959接下来研究含有未知参数的本例中，''(4cos2)015yqxyq，''(0)1.(0)0,()0yyy微分方程与参数λ的数值有关。一般而言，对于任意的λ值，该问题无解，但对于特殊的λ值（特征值），它存在一个解，这也称为微分方程的特征值问题。对00.10.20.30.40.50.60.70.80.9111.21.41.61.822.22.42.62.83AnalyticalSolutionNumericalSolution于此问题，可在bvpinit中提供参数的猜测值，然后重复求解BVP得到所需的参数，返回参数为sol.parameters。functionb_bvp2lamb=15;solinit=bvpinit(linspace(0,pi,10),[1;1],lmb);opts=bvpset('Stats','on');sol=bvp4c(@ODEfun,@BCfun,solinit,opts);lamb=sol.parametersx=[0:pi/60:pi];y=deval(sol,x);plot(x,y(1,:),'b-',sol.x,sol.y(1,:),'ro')legend(''解曲线','初始网格点解')%------------------------------------------------------------------functiondydx=ODEfun(x,y,lamb)q=5;dydx=[y(2);-(lamb-4*q*cos(2*x))*y(1)];%------------------------------------------------------------------functionbc=BCfun(ya,yb,lmb)bc=[ya(1)-1;ya(2);yb(2)];于是在commandwindow中运行b_bvp2可以得到：Thesolutionwasobtainedonameshof53points.Themaximumresidualis6.960e-04.Therewere1701callstotheODEfunction.Therewere64callstotheBCfunction.lmb=21.1051改变参数q=9结果如下：lmb=26.9888结果分析与讨论：从上面的结果中可以看出，对于第一个微分方程组使用bvp4c00.511.522.533.5-0.8-0.6-0.4-0.200.20.40.60.811.2解曲线初始网格点解00.511.522.533.5-1.5-1-0.500.511.52解曲线初始网格点解函数的求解，比较简单，直接根据微分方程的形式编写函数文件，要注意bvp4c的调用规则，同时要注意初值边界条件。把待解的问题转化为标准边值问题，因为边值问题可以多解，所以需要为期望解指定一个初始猜测解。该猜测解网（Mesh）包括区间[a,b]内的一组网点（Meshpoints）和网点上的解S(x)根据原微分方程构造残差函数，利用原微分方程和边界条件，借助迭代不断产生新的S(x)，使残差不断减小，从而获得满足精度要求的解。对于第二种情况，微分方程与参数λ的数值有关。一般而言，对于任意的λ值，该问题无解，但对于特殊的λ值（特征值），它存在一个解，这也称为微分方程的特征值问题。对于此问题，可在bvpinit中提供参数的猜测值，然后重复求解BVP得到所需的参数，返回参数为sol.parameters。这是需要特别注意的地方。