数学建模社区-数学中国

标题: 线性微分方程边值问题打靶算法Matlab程序 [打印本页]
作者: 建不了的模。    时间: 2014-10-6 10:35
标题: 线性微分方程边值问题打靶算法Matlab程序
注意该算法只能完成二阶常微分方程双边值问题求解，至于其他形式的边值问题必须先转换到二阶形式

对于下面的二阶常微分方程

20090319_ee1f3a88f0e051ac1a7bZlH10CkxnFuy.gif

file:///http://attach.matlabsky.com/data/attachment/forum/month_0903/20090319_ee1f3a88f0e051ac1a7bzlh10ckxnfuy.gif
利用上面方程的线性结果和两个特殊的初值问题，我们可以构造两个等效的常微分初值方程初值问题，对于初值问题我们就可以直接使用ode**计算器或者龙哥库塔算法求解了。构造的两个初值问题为

20090319_aef647afd9e5f3531724bbuILIol73n3.gif

此时我们可以构造原微分方程的解为
20090319_317956a9be701da19363jqIYCxERXeok.gif
有上面的基础那我们很容编写直接边值问题的线性打靶算法，下面我写了一个供大家参考(参见附件)，函数内部调用了自己编写的rk4()函数
%线性方程边值问题的打靶算法
%对微分方程y''=p*y'+q*y+r,a<t<b,u(a)=alpha,u(b)=beta,其中y,p,q和r都是x的函数
%构造成两个等效的初值问题
%u''=p*u'+q*u+r,u(a)=alpha,u'(a)=0
%v''=p*v'+q*v,v(a)=beta,v'(a)=1
%使用Matlab语言描述微分方程u和v，注意我们先必须化为一阶微分方程组
%f1=@(t,u)[u(2);p*u(2)+q*u(1)+r];
%f2=@(t,v)[v(2);p*v(2)+q*v(1)];
%则原微分方程的解为
%x=u+c*v=u+(beta-u(b))[img]file:///C:\DOCUME~1\ADMINI~1\LOCALS~1\Temp\[LC3U)F{0XCAB)LKNIT0K@G.gif[/img](b)*v
%
%输入参数
% f1 描述方程u的Matlab函数,只能是句柄(@()和M函数)，或者inline函数
% f2 描述方程v的Matlab函数
% a,b 微分计算区间[a,b]
% alpha,beta 上下边界
% maxiter 最大迭代次数
% 输出参数
% [t,x] 坐标点t对应的微分值x
%
我们举一个例子吧，对于下面的二阶边值问题

20090319_26b5bc282ca5d1a96d2fNt6OK1vPj0kU.gif

f1=@(t,u)[u(2);2*t/(1+t.^2)*u(2)-2/(1+t.^2)*u(1)+1];%这里必须使用“;”隔开
f2=@(t,u)[u(2);2*t/(1+t.^2)*u(2)-2/(1+t.^2)*u(1)];
a=0;b=4;alpha=1.25;beta=-0.95;maxiter=1000;
[t,x]=lineshoot(f1,f2,a,b,alpha,beta,maxiter);
plot(t,x)
复制代码

20090319_91b2a46deb5b9a7a0fd3oIaDVwXjqUhX.jpg.thumb.jpg
附件中的代码使用到了自己编写的四阶龙阁库塔算法rk4()函数，但是假如说对Matlab中的ode**函数比较熟悉的网友，可以直接调用ode45。

下面给出参考代码，但是程序是将原来的边值问题转换为如下四个初值问题：

20090319_6283cb31a95044720aa1uob2NLr6kQFz.gif
lineshoot.rar (1.68 KB, 下载次数: 0)
作者: 深V礼    时间: 2014-10-6 13:34
不错的书籍，下载来看看
作者: ゞ_轻描丶幸福的    时间: 2014-10-6 13:35
好资料   下载看看
作者: mingtingqing    时间: 2014-10-6 13:36
太棒了
作者: fanfan1993    时间: 2015-2-3 13:30
谢谢~~~~~~~
作者: jiay    时间: 2015-2-3 14:23
good............




欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5