前置知识-边值问题、打靶法、bvp 系列函数的用法

1.2 边值问题

微分方程边值问题（Boundary Value Problem，简称BVP）是微分方程求解中的一个重要问题。与初值问题（Initial Value Problem，简称IVP）不同，BVP是在某个区间内寻求微分方程解的特定边界条件下的解。

在实际问题中，许多微分方程的解必须满足一些特定的边界条件，例如固定的温度或者压力等。这时就需要使用BVP求解微分方程的解。BVP的解通常比IVP更加复杂，因为它需要同时满足微分方程和一组特定的边界条件。

BVP的解决方法通常需要使用一些数值方法，例如有限差分法、有限元法和谱方法等。这些方法的主要思想是将微分方程转化为一个线性代数方程组，然后使用矩阵运算求解。这些方法的优点是可以处理一般的非线性和非均匀微分方程，但是它们需要大量的计算和数值处理，并且需要选择合适的数值方法和参数。

BVP的应用范围非常广泛，例如在工程学中可以用于处理热传导、流体力学、结构力学等问题，在物理学中可以用于处理量子力学、电动力学等问题。因此，了解和掌握BVP求解方法对于理解和应用微分方程有着重要的意义。

本小节讨论典型的二阶常微分方程的边值问题:
$$\frac{{\mathrm{d}}^{2}u}{\mathrm{d}{x}^2}=f\left(x,u,\frac{\mathrm{d}u}{\mathrm{d}x}\right),a<x<b\text{\hspace{1em}(1-65)}$$
其边界条件为 $({\alpha }_1、{\beta }_1、{\gamma }_1、{\alpha }_2、{\beta }_2、{\gamma }_2$ 均为常数 $)$ :
$$\{\begin{array}{l}{\alpha }_{1}u(a)+{\beta }_1{u}^{\prime }(a)={\gamma }_1\\ {\alpha }_{2}u(b)+{\beta }_2{u}^{\prime }(b)={\gamma }_2\end{array}\text{\hspace{1em}(1-66)}$$
上式具有很强的一般性, 这里先从一种简单的边界条件开始分析 ( $\alpha 、\beta$ 为常数):
$$u(a)=\alpha ,u(b)=\beta \text{\hspace{1em}(1-67)}$$

1.2.1 打靶法.

对于二阶常微分方程（1-65）及边界条件（1-67), 利用代换 $u_1=u、u_2=u^{\prime }$, 可将其写为等价的方程组形式:
$$\{\begin{array}{ll}{u}_1^{\prime }=u_2,& u_1(a)=\alpha \\ u_2^{\prime }=f\left(x,u_1,u_2\right),& u_1(b)=\beta \end{array}\text{\hspace{1em}(1-68)}$$
它与初值问题的不同之处在于, $u_2(a)$ 的值是末知的, 已知条件中取而代之的是 $u_1(b)$ 的值, 否则就可以直接通过 odesolver 数值求解了。假如可以确定一个值 $m$, 使初值问题 (1-69) 的解 $u_1(x)$ 满足 $\mid u_1(b)-\beta \mid <\epsilon$, 其中 $\epsilon$ 为可允许的误差, 那么边值问题 (1-68) 和初值问题（1-69) 的解就是近似相等的, 也就相当于把边值问题转化成了初值问题。
$$\{\begin{array}{ll}{u}_1^{\prime }=u_2,& u_1(a)=\alpha \\ u_2^{\prime }=f\left(x,u_1,u_2\right),& u_2(a)=m\end{array}\text{\hspace{1em}(1-69)}$$
打靶法就是用来确定 $m$ 值的, 大致思路如下:

(1) 随意选取 $m$ 的起始值。

(2) 求出初值问题 (1-69) 的解 $u_1(x)$ 。

(3）若不符合 $\mid u_1(b)-\beta \mid <\epsilon$, 就修正 $m$ 的值。

(4) 重复前面两项, 直到符合 $\mid u_1(b)-\beta \mid <\epsilon$, 此时的 $u_1(x)$ 就是数值解。

修正 $m$ 取值的方法可以是线性插值法、牛顿法等, 这里使用比较简单的线性插值法。 设 $m=m_n$ 时初值问题（1-69) 的解 $u_1(x)$ 在 $x=b$ 处有 $u_1(b)={\beta }_n$, 其中 $n=1,2,\dots$ 。若假设 $m_n$ 与 ${\beta }_n$ 是线性关系, 为了让 $m$ 取最新修正值 $m_{n+2}$ 时有 ${\beta }_{n+2}=\beta$, 则 $m_{n+2}$ 需要满足 $\left(m_{n+2}-m_n\right)/\left(\beta -{\beta }_n\right)=\left(m_{n+1}-m_n\right)/\left({\beta }_{n+1}-{\beta }_n\right)$, 即:
$$m_{n+2}=m_n+\frac{m_{n+1}-m_n}{{\beta }_{n+1}-{\beta }_n}\left(\beta -{\beta }_n\right)\text{\hspace{1em}(1-70)}$$
其中的 $m_n、{\beta }_n、m_{n+1}、{\beta }_{n+1}$ 为已知, 这就是说如果使用线性插值法修正 $m$ 的取值, 在打靶法开始时需要先随意选取 2 个 $m$ 的起始值 $m_1、m_2$, 计算出相应的 ${\beta }_1、{\beta }_2$, 才可以根据上式得到修正值 $m_3$, 进而得到 ${\beta }_3$, 然后得到修正值 $m_4\cdots \cdots$ 直到找到满足 $\mid {\beta }_n-\beta \mid <\epsilon$ 的 $m$ 值。此时初值问题 (1-69) 的解即近似为边值问题（1-68）的解, 而初值问题是容易数值求解的。

虽然 $m$ 取 $m_1、m_2$ 时, 最终结果在 $x=b$ 处的取值与应有的边界条件相差甚远, 但通过一次一次地修正, 逐渐缩小距离目标的差距, 最终总会得到在一定误差内满足 $x=b$ 处边界条件的解。

对于更复杂的边界条件 (1-66), 照样可以用打靶法解决。尽管 $u(a)$ 和 $u^{\prime }(a)$ 的值都是末知的, 但只要猜测它们其中一个的取值, 就可以通过关系 ${\alpha }_{1}u(a)+{\beta }_1{u}^{\prime }(a)={\gamma }_1$ 得到另外一 个的值。有了 $u(a)$ 和 $u^{\prime }(a)$ 的值, 再通过 odesolver 求解式 (1-65) 得到 $u(b)$ 和 $u^{\prime }(b)$ 的值（当然需要变量代换 $u_1=u、u_2=u^{\prime }$ ), 然后使用 $\mid {\alpha }_{2}u(b)+{\beta }_2{u}^{\prime }(b)-{\gamma }_2\mid <\epsilon$ 判断猜测的初始值是否符合 要求, 如失败就修正猜测值后继续尝试, 直到符合要求。这样一来, 同样可以将边界条件 （1-66）下的边值问题当做初值问题处理。

下面举一个例子, 考虑如下二阶常微分方程:
$$u^{\prime \prime }+|u|=0,0<x<4\text{\hspace{1em}(1-71)}$$
其边界条件为:
$$u(0)=0,u(4)=-2\text{\hspace{1em}(1-72)}$$
做代换 $u_1=u、u_2=u^{\prime }$, 写为一阶常微分方程组形式:
$$\{\begin{array}{ll}{u}_1^{\prime }=u_2,& u_1(0)=0\\ u_2^{\prime }=-\left|u_1\right|,& u_1(4)=-2\end{array}\text{\hspace{1em}(1-73)}$$
将此问题转化为初值问题, 就需要找到一个 $m=m_n$, 使式 (1-74) 的解 $u_1(x)$ 在 $x=4$ 处 的值 $u_1(4)={\beta }_n$ 满足 $\mid {\beta }_n+2\mid <10^{-6}$ 。
$$\{\begin{array}{ll}{u}_1^{\prime }=u_2,& u_1(0)=0\\ u_2^{\prime }=-\left|u_1\right|,& u_2(0)=m\end{array}\text{\hspace{1em}(1-74)}$$
根据打靶法思想, 确定 $m$ 取值的流程如下所示

相应代码如下:

主程序代码如下

clear ; close all;
n=0; m(1)=1; alpha=0; beta0=-2;
while n==0||abs(beta(n)-beta0)>=1e-6
    n=n+1;
    %在猜测的初始条件下求解初值问题
    [x,usol]=ode45('shoot',[0 4],[alpha m(n)]);
    beta(n)=usol(end,1);
    if n==1
        m(2)=m(1)-0.1;
    else
        %线性插值法修正m的取值
        m(n+1)=m(n-1)+(m(n)-m(n-1))*(beta0-beta(n-1))/(beta(n)-beta(n-1));
    end
end
plot(x,usol(:,1),'LineWidth',1.5), xlabel x, ylabel u

文件 shoot.m 代码如下

function du=shoot(x,u)
du=[u(2); -abs(u(1))];

代码中使用的 $m$ 的前两个猜测值为 1 和 $0.9$, 程序仅经过了 5 次循环, 就满足条件 $\mid {\beta }_5+2\mid <10^{-6}$, 得到结果如下图所示。如选取 $m$ 的前两个猜测值为-1 和-1.1, 还会得到另 外一个截然不同的解。因为边值问题的解往往不是唯一的, 所求得的解与猜测值的选取有关, 这是边值问题与初值问题的一个很大不同。

查看程序运行后 $m_n$ 和 ${\beta }_n$ 的前 5 个值, 可发现线性插值法的修正效率还比较令人满意。

>> m(1:5)

ans =

    1.0000    0.9000    2.0480    2.0640    2.0641

>> beta(1:5)

ans =

   -0.9699   -0.8716   -1.9845   -1.9999   -2.0000

1.2.2 bvp 系列函数的用法

与专门针对初值问题的 ode 系列函数类似, Matlab 中还有专门用于解决边值问题的 bvp 系列函数, 其调用语法为:

sol = bvpsolver(odefun, bcfun,solinit,options)

下面对每一项进行说明:

(1) “bvpsolver” 是 bvp 系列函数的名称, 可以是 bvp4c 或 bvp5c, 后者比前者的精度高一些。

(2) “odefun" 是包含常微分方程组的函数句柄, 这与 ode 系列函数中的 odefun 是一样的。

(3) “bcfun" 是包含边界条件的函数句柄。

(4) “solinit”是包含初始网格结点、初始猜测值的结构体, 由 bvpinit 函数生成。

(5) “options”用于改变 bvpsolver 的默认设置, 由 bvpset 函数生成, 为可选项。

(6) “sol"为包含计算结果的结构体。

如不需要用到 options, 则 bvp 系列函数的调用形式为:

sol = bvpsolver(odefun,bcfun,solinit)

用于生成结构体 solinit 的 bvpinit 函数的调用语法为:

solinit=bvpinit(x,uinit,params)

其中:

(1) “ x ”为计算区间内的初始网格结点, 起止点即为边界条件所在的两处位置, 这些结点需要按大小单调排列。

(2) “uinit" 是解的初始猜测值, 也可以是每个结点处的解的猜测值。

(3) “params" 是边值问题中的末知参数 (如果有的话) 的猜测值, 为可选项。

下面举例说明 bvp4c 的用法, 边值问题如下:
$$\{\begin{array}{l}{u}^{\prime \prime }=-u+2(x-1){\mathrm{e}}^{-x},0<x<\pi \\ u(0)=1,u(\pi )=\pi /{\mathrm{e}}^{\pi }-1\end{array}\text{\hspace{1em}(1-75)}$$
先利用代换 $u_1=u、u_2=u^{\prime }$, 化为一阶常微分方程组形式:
$$\{\begin{array}{ll}{u}_1^{\prime }=u_2,& u_1(0)=1\\ u_2^{\prime }=-u_1+2(x-1){\mathrm{e}}^{-x},& u_1(\pi )=\pi /{\mathrm{e}}^{\pi }-1\end{array}\text{\hspace{1em}(1-76)}$$
将此常微分方程组和边界条件分别写进 shoot $2.\mathrm{m}$ 和 bc.m 中, 用 bvp4c 求解, 代码如下:

主程序代码如下

clear; close all;
x=linspace(0,pi,20);
solinit=bvpinit(x,[2 2]);
%bvp4c解边值问题
sol=bvp4c(@shoot2,@bc,solinit);
%从结构体sol获取指定位置的数值结果
u=deval(sol,x);
%解析解
u_exact=cos(x)+x.*exp(-x);
%画图
plot(x,u(1,:),'+k',x,u_exact,'MarkerSize',10,'LineWidth',1.5)
set(gca,'Fontsize',16), xlabel x, ylabel u
legend('bvp4c','解析解'), axis([0 pi -1 1.5])

文件 bc.m 代码如下

function res=bc(ua,ub)
%边界条件boundary condition
res=[ua(1)-1; ub(1)+1-pi/exp(pi)];

这里重点说明一下包含边界条件的文件 bc.m, 其中定义了一个同名函数 bc, 参数 ua 和 ub 分别代表左边界和右边界的情况。其中, ua(1) 和 ua(2) 为函数在左边界的取值和一阶导数, ub(1) 和 ub(2) 为函数在右边界的取值和一阶导数。代码中的 ua(1)-1 就是 $ua(1)-1=0$ 的简写, 代表边界条件 $u(0)=1$, 同理, ub(1)+1-pi/exp(pi)表示 $u(\pi )=\pi /{\mathrm{e}}^{\pi }-1$ 。

代码中的初始猜测值为 $u_1=2、u_2=2$, 假如边值问题的解是唯一的, 初始猜测值的选取是不会影响结果的。下图同时显示了 bvp4c 得到的数值解和式 (1-75) 的解析解 $u=\cos (x)+x{\mathrm{e}}^{-x}$, 二者相当吻合。值得一提的是, 如果在 bvpinit 函数中选取的初始网格结点数量偏少的话, 会造成数值解和解析解存在一定的差异, 在这一点上, bvp4c 比打靶法更难使用。