打靶法求解两点边值问题简介及Matlab程序

weixin

　　常微分方程是我们在解决实际问题中经常遇到的数学模型，在具体求解常微分方程时，必须附加某种定解条件。定解条件通常有两种，一种是初始条件，另一种是边界条件，与边界条件相应的定解问题称为边值问题。本文主要对求解两点边值问题的打靶法进行简要的介绍。

　　对于工程中常见的方程：

　　当f关于y和y'是线性时，上式为线性两点边值问题，可表示为：

　　打靶法的基本原理是将两点边值问题转化为下列形式的初值问题：

　　这里的sk为y在a处的斜率。令z=y'，上述二阶方程可降为一阶方程组：

　　因此，边值问题变成求合适的sk，使上述方程组初值问题的解满足原边值问题的右端边界条件y(b)=β，从而得到边值问题的解。这样，把一个两点边值问题的数值解问题转化为一阶方程组初值问题的数值解问题。方程组初值问题的所有数值方法在这里都可以使用。问题的关键是如何去找合适的初始斜率的试探值sk。

　　对给定的sk，设初值问题的解为y(x,sk)，它是sk的隐函数。假设y(x,sk)随sk是连续变化的，记为y(b,s)-β=0，于是我们要找的sk就是方程的根。可以用迭代法求上述方程的根。比如用割线法有：

　　这样，可以按下面简单的计算过程进行求解。先给定两个初始斜率s0、s1，分别作为初值问题的初始条件。用一阶方程组的数值方法求解它们，分别得到区间右端点的函数的计算值y(b,s0)和y(b,s1)。如果：

　　则以y(b,s0)或y(b,s1)作为两点边值问题的解。否则用割线法求s2，同理得到y(b,s2)，再判断它是否满足精度要求。

　　如此重复，直到某个sk满足：

　　此时得到的y(xi)和y'i=z(xi)就是边值问题的解函数值和它的一阶导数值。上述方程好比打靶，sk作为斜率为子弹的发射，y(b)为靶心，故称为打靶法。

　　值得指出的是，对于线性边值问题，一个简单又实用的方法是用解析的思想，将它转化为两个初值问题：

　　求得这两个初值问题的解y1(x)和y2(x)，若y2(b)≠0，容易验证线性两点边值问题的解为：

　　以下为新浪博客博主xianfa110分享的打靶法求解两点边值问题的Matlab程序，仅供参考。

　　function varargout = shooting_two_point_boundary(varargin)
　　% =====================
　　% 函数名：shooting_two_point_boundary.m
　　% 基于打靶法计算两点边值问题
　　% 仅针对二阶微分方程
　　% author: xianfa110.
　　% blog.sina.com.cn/xianfa110
　　% 函数形式：
　　% [result,err,z0] = shooting_two_point_boundary(@fun,[y_0,y_end],[x_0,x_1],h);
　　% 输入：
　　% fun = 函数名;
　　% y_0 = 函数初值;
　　% y_end = 函数终值;
　　% x_0 = 自变量初值;
　　% x_end = 自变量终值;
　　% h = 积分步长;
　　% 输出：
　　% result = [x,y];
　　% err = 误差;
　　% z0 = y'初值;
　　% =====================
　　% 函数fun：4y''+yy' = 2x^3 +16 ; 2<= x <=3
　　% 写法：
　　% function f = fun(y,x)
　　% dy = y(2);
　　% dz = (2*x^3+16-y(1)*y(2))/4;
　　% f = [dy,dz];
　　% =====================
　　% 注意:y(1) = y，y(2) = y'。
　　% =====================
　　F = varargin{1};
　　y_0 = varargin{2}(1);
　　y_end = varargin{2}(2);
　　x_0 = varargin{3}(1);
　　x_1 = varargin{3}(2);
　　ts = varargin{4};
　　t0 = x_0-0.5;
　　flg = 0;
　　kesi = 1e-6;
　　y0 = rkkt(F,[y_0,t0],x_0,x_1,ts);
　　n = length(y0(:,1));
　　if abs(y0(n,1)-y_end)<=kesi
　　   flg=1;
　　else
　　   t1=t0+1;
　　   y1=rkkt(F,[y_0,t1],x_0,x_1,ts);
　　   if abs(y1(n,1)-y_end)<=kesi
　　      flg=1;
　　   end
　　end
　　if flg ~= 1
　　   while abs(y1(n,1)-y_end) > kesi
　　   % ==插值法求解非线性方程==
　　      t2 = t1-(y1(n,1)-y_end)*(t1-t0)/(y1(n,1)-y0(n,1));
　　      y2 = rkkt(F,[y_0,t2],x_0,x_1,ts);
　　      t0=t1;
　　      t1=t2;
　　      y0=y1;
　　      y1=y2;
　　   end
　　end
　　x = x_0:ts:x_1;
　　out = [x',y1(:,1)];
　　varargout{1} = out;
　　varargout{2} = abs(y1(n,1)-y_end);
　　varargout{3} = t1;

　　本文根据网络资料整理而成。