实验八_常微分方程边值问题.doc

实验八—边值问题的求解 在2011年全国大学生数学建模邀请赛的A题的数学模型可归结为如下常微分方程组的边值问题: (1) 其中，， 是参数. 边值条件是 (2) 从物理意义上考虑，可分成以下3种情况来讨论: 1当梁不与地面接触时，要去掉(2)中最后一个边值条件. 并且c=1，p=1，四个边值条件就惟一确定了解, 记，. .2当梁仅与地面是线接触时，去掉(2)中最后一个边值条件. 并且c=1，这时存在一个临界值，，使得当时，四个边值条件惟一确定了解，当从1单调减到时，从单调减少到零. 从单调减少到某个.故当时，存在惟一的，使得. 从而解决了情况2. 而当p从进一步减少时，将会小于零，不再符合与水平地面接触的物理条件. 而应该考虑以下“面接触”的情况. 3当梁与地面是面接触时，要考虑(2)中全部五个边值条件，当确定时，五个边值条件就惟一确定了c1, 及[0，c]上方程的解. 当p从单调减少到0时，从单调减少到0. 故当时，存在惟一的，使得. 从而解决了情况3. c是个与有关的，故这个问题是个活动边界问题，它比情况1，2的普通的边值问题更难求解. 而当参数a=24，h=1/4时就属于第3种情况. 普通的边值问题求解的基本方法是打靶法. 对于第1种情况，打靶法就是考虑初值，选择不同的，求解常微分方程组(1), 使得. 对于第2种情况，就是考虑初值，选择不同的及p，求解常微分方程组(1), 使得. ， 对于第3种情况，就是考虑初值，选择不同的及p，c，求解常微分方程组(1), 使得. ，. 求解方法1：以上三种情况的求解皆可以通过ode45求解，得到一个MATLAB函数，输入是初值及参数中的未知量，输出是终值减去条件值，然后数值求解这个MATLAB函数零点即可得到初值及参数中的未知量，再数值求解常微分方程组的初值问题的解即可. 然而求根时初始值的选取是很关键的. 对于第1种情况，就是求一个一元函数的根，可以通过画图得到根的近似值作为初始值. 对于第2种情况，就是求一个二元函数的根的问题，考虑到根应当是h的连续函数，故可以通过逐渐减小h的值来得到根的近似初始值. 具体来说，就是如对于某个值，求出了函数的根，把这个根作为对于较小的h时函数的根的近似，照此一直进行，直到求出实际的时函数的根. 同理用这种方法可以求解第3种情况. 在求解之前，必须要求出才能确定是哪种情况. 对于第1种情况还可以直接利用MATLAB求边值问题的函数bvp4c来求解，bvp4c的完全的用法为 sol = bvp4c(@Odefun, @Bcfun,Solinit,Options,p1,p2...) 其中Odefun 是微分方程函数，形式为 dydx = odefun(x,y,parameters,p1,p2,...) 其中x是自变量, y is 未知的列向量. parameters 是未知的参数向量, p1,p2,... 是已知的参数. dydx 是列向量. Bcfun 是自编的边值条件的残差的函数，对于二点边值问题，函数Bcfun的格式(例如) function res = Bcfun(ya,yb,parameters,p1,p2,...) ya=[y(1), y(2),y(3)]; yb=y(4); 其中ya是左边值条件中的残差，yb是右边值条件的残差 Solinit 是你对解的初始猜测. 是通过函数bvpinit 得到的. Solinit=bvpinit(xinit,yinit，params); 其中xinit是自变量的点列，yinit是猜测的对应的因变量的点列，yinit(:,i) 是在x(i)处的解的列向量 params是猜测的parameters的值. 猜测的值在边界上最好是满足边值条件. 当以上parameters,p1,p2之中没有的可以缺省不写. Options是通过函数bvpset设置精度等的选项. 如 Options=bvpset(RelTol,5.04968e-013,AbsTol,6.33975e-013); 而默认的RelTol=1e-3, AbsTol=1e-6 对于第一小问，以下整个程序(可以保存为一个m文件bvp.m)的功能为给定一个a0的值，用bvp4c求解边值问题,输出右端的y值h，并画出梁的图像 function h=bvp(a) Yinit=@(x) Yini(x,a); %初始猜测的解 Solinit = bvpinit(linspace(0,1,10),Yinit); %初始猜测的解的数据点 Options=bvpset(RelTol,5.04968e-013,AbsTol,6.33975e-013); %精度设置 sol = bvp4c