本页对应的英文页面已更新,但尚未翻译。 若要查看最新内容,请点击此处访问英文页面。

bvp4c

求常微分方程的边界值问题的解

语法

sol = bvp4c(odefun,bcfun,solinit)
sol = bvp4c(odefun,bcfun,solinit,options)
solinit = bvpinit(x, yinit, params)

参数

odefun

计算微分方程 f(x,y) 的函数句柄。其形式可以为

dydx = odefun(x,y)
dydx = odefun(x,y,parameters)

对于标量 x 和列向量 yodefun(x,y) 必须返回列向量 dydx,表示 f(x,y)。parameters 是一个未知参数向量。

bcfun

计算边界条件中的残差的函数句柄。对于 bc(y(a),y(b)) 形式的两点边界值条件,bcfun 的形式可以为

res = bcfun(ya,yb)
res = bcfun(ya,yb,parameters)

其中 yayb 是与 y(a)y(b) 对应的列向量。parameters 是未知参数的向量。输出 res 是列向量。

有关适用于多点边界值问题的 bcfun 的说明,请参阅多点边界值问题

solinit

包含初始估计解的结构体。应使用函数 bvpinit 创建 solinitsolinit 具有以下字段。

 

x

初始网格的排序节点。边界条件应用在a = solinit.x(1) 且 b = solinit.x(end) 处。

 

y

初始估计解,使得 solinit.y(:,i) 为在 solinit.x(i) 节点处的初始估计解。

 

parameters

可选。提供未知参数的初始估计值的向量。

 

结构体可以随意命名,但是字段必须命名为 xyparameters。可以使用辅助函数 bvpinit 构建 solinit。有关详细信息,请参阅 bvpinit

options

可选积分参数。使用 bvpset 函数创建的结构体。有关详细信息,请参阅 bvpset

说明

sol = bvp4c(odefun,bcfun,solinit) 计算以下形式的常微分方程组的积分

y′ = f(x,y)

(在取决于以下两点边界值条件的区间 [a,b] 中 bc(y(a),y(b)) = 0)。

odefunbcfun 均为函数句柄。有关详细信息,请参阅创建函数句柄

参数化函数 介绍如何在必要情况下向函数 odefun 和边界条件函数 bcfun 提供其他参数。

bvp4c 还可以求多点边界值问题的解。请参阅 多点边界值问题。使用函数 bvpinit 可以指定边界点,这些边界点存储在输入参数 solinit 中。有关详细信息,请参阅 bvpinit 的参考页。

bvp4c 求解器还可以查找以下形式的问题的未知参数 p

y′ = f(x,y, p)

0 = bc(y(a),y(b),p)

其中 p 与 parameters 对应。应在 solinit.parameters 中为 bvp4c 提供任意未知参数的初始估计值。bvp4c 求解器会在 sol.parameters 中返回未知参数的最终值。

bvp4c 会生成在 [a,b] 上的连续解,在其中该解具有连续一阶导数。使用函数 devalbvp4c 的输出 sol 来计算区间 [a,b] 中的特定点 xint 的解。

sxint = deval(sol,xint)

bvp4c 返回的结构体 sol 包含下列字段:

sol.x

bvp4c 选择的网格

sol.y

sol.x 网格点处的 y(x) 近似值

sol.yp

sol.x 网格点处的 y′(x) 近似值

sol.parameters

bvp4c 针对未知参数返回的值(如果有)

sol.solver

'bvp4c'

sol.stats

计算成本统计信息(当设置 stats 选项和 bvpset 时,也会显示此信息)。

结构体 sol 可以随意命名,bvp4c 创建字段 xyypparameterssolver

sol = bvp4c(odefun,bcfun,solinit,options) 的解算方法与上述方法相同,但将默认积分属性替换为 options(使用 bvpset 函数创建的结构体)中的值。有关详细信息,请参阅 bvpset

solinit = bvpinit(x, yinit, params) 使用未知参数的估计值的向量 params 得出初始估计值 solinit

奇异边界值问题

bvp4c 可对一类奇异边界值问题求解,其中包括具有未知参数 p 的以下形式的问题:

y′ = S · y/x + F(x,y,p)

0 = bc(y(0),y(b),p)

需要区间为 [0, b] 并且 b > 0。当计算偏微分方程 (PDE) 因柱状或球面对称性而产生的 ODE 的平滑解时,通常会出现此类问题。对于奇异问题,应将(常量)矩阵 S 指定为 bvpset'SingularTerm' 选项的值,并且 odefun 仅计算 f(x,y,p)。边界条件必须与必要条件 S·y(0) = 0 一致,并且初始估计值应满足此必要条件。

多点边界值问题

bvp4c 可以求解多点边界值问题(其中 a = a0 < a1 < a2 < ...< an = b 是区间 [a,b] 中的边界点)。点 a1,a2,...,an−1 表示将 [a,b] 分为多个区域的界点。bvp4c 使用从 1 开始的索引按从左到右的顺序(从 a 到 b)枚举区域。在区域 k [ak−1,ak] 中,bvp4c 按如下所示计算导数

yp = odefun(x,y,k) 

在以下边界条件函数中

bcfun(yleft,yright) 

yleft(:,k) 是 [ak−1,ak] 左边界的解。同样,yright(:,k) 是区域 k 右边界的解。具体来说是

yleft(:,1) = y(a)

yright(:,end) = y(b)

当使用以下公式创建初始估计值时,

solinit = bvpinit(xinit,yinit),

请在每个接口点的 xinit 中使用双精度输入。有关详细信息,请参阅 bvpinit 的参考页。

如果 yinit 为函数,bvpinit 调用 y = yinit(x,k) 以获取区域 k 中的 x 处的初始估计解。在 bpv4c 返回的解结构体 sol 中,sol.x 具有每个接口点的双精度输入。sol.y 的对应列分别包含界点的左解和右解。

要查看对三点边界值问题求解的示例,请在 MATLAB® 命令提示符下键入 threebvp

注意

除两个求解器之间的误差容限含义以外,bvp5c 函数的使用方法与 bvp4c 完全类似。如果 S(x) 近似于解 y(x),则 bvp4c 控制残差 |S′(x) – f(x,S(x))|。这会间接控制真误差 |y(x) – S(x)|。bvp5c 直接控制真误差。

示例

使用初始估计值指示所需解

边界值问题可能有多个解。初始估计值的一个目的是指示您需要多个解中的哪个解。

二阶微分方程

实际上有两种满足边界条件的解

在使用 bvp4c 解决该问题之前,您需要将该微分方程重写为包含两个一阶 ODE 的方程组,

其中 并且 。此方程组具有所要求的形式

函数 f 和边界条件 bc 在 MATLAB 中编码为函数 twoodetwobc

function dydx = twoode(x,y)
dydx = [ y(2); -abs(y(1)) ];
function res = twobc(ya,yb)
res = [ ya(1); yb(1) + 2 ];

构造估计值结构体,该结构体包含 [0,4] 区间内五个等距点的初始网格以及以下估计常量值

solinit = bvpinit(linspace(0,4,5),[1 0]);

使用 bvp4c 解决该问题。

sol = bvp4c(@twoode,@twobc,solinit);

以 100 个等间距点计算数值解并对 y(x) 绘图。

x = linspace(0,4);
y = deval(sol,x);
plot(x,y(1,:))

要获得此问题的其他解,请使用初始估计值

solinit = bvpinit(linspace(0,4,5),[-1 0]);
sol = bvp4c(@twoode,@twobc,solinit);
x = linspace(0,4);
y = deval(sol,x);
plot(x,y(1,:))

计算 Mathieu 方程的第四个特征值

此边界值问题涉及未知参数。此任务的目的是计算以下 Mathieu 方程的第四个 (q = 5) 特征值 λ

y” + (λ – 2q cos2x)y = 0。

因为存在未知参数 λ,所以此二阶微分方程受限于三个边界条件:

y′(0) = 0

y′(π) = 0

y(0) = 1

使用局部函数将 bvp4c 所需的所有函数放置在一个文件中非常方便。

function mat4bvp

lambda = 15;
solinit = bvpinit(linspace(0,pi,10),@mat4init,lambda);
sol = bvp4c(@mat4ode,@mat4bc,solinit);

fprintf('The fourth eigenvalue is approximately %7.3f.\n',...
        sol.parameters)

xint = linspace(0,pi);
Sxint = deval(sol,xint);
plot(xint,Sxint(1,:))
axis([0 pi -1 1.1])
title('Eigenfunction of Mathieu''s equation.') 
xlabel('x')
ylabel('solution y')
% ------------------------------------------------------------
function dydx = mat4ode(x,y,lambda)
q = 5;
dydx = [  y(2)
         -(lambda - 2*q*cos(2*x))*y(1) ];
% ------------------------------------------------------------
function res = mat4bc(ya,yb,lambda)
res = [  ya(2) 
         yb(2) 
        ya(1)-1 ];
% ------------------------------------------------------------
function yinit = mat4init(x)
yinit = [  cos(4*x)
          -4*sin(4*x) ];

微分方程(转换为一阶方程组)和边界条件分别编码为局部函数 mat4odemat4bc。由于存在未知参数,因此这些函数必须接受三个输入参数,即使未使用某些参数也是如此。

估计值结构体 solinit 使用 bvpinit 构造而成。解的初始估计值以函数 mat4init 形式提供。我们选择 y = cos 4x 是因为它满足边界条件,并且具有正确的定性行为(正确的符号更改数)。在调用 bvpinit 时,第三个参数 (lambda = 15) 会提供未知参数 λ 的初始估计值。

在使用 bvp4c 对该问题求解之后,字段 sol.parameters 返回值 λ = 17.097,并且绘图显示与此特征值关联的特征函数。

算法

bvp4c 是一个有限差分代码,此代码实现 3 阶段 Lobatto IIIa 公式。这是配置公式,并且配置多项式会提供在整个 [a,b] 中具有一致四阶精度的 C1 连续解。网格选择和误差控制均基于连续解的残差。

参考

[1] Shampine, L.F., M.W. Reichelt, and J. Kierzenka, “Solving Boundary Value Problems for Ordinary Differential Equations in MATLAB with bvp4c

在 R2006a 之前推出