# bvp4c
求解边界值问题 - 四阶方法
函数库: TyMath
# 语法
sol = bvp4c(odefun,bcfun,solinit)
sol = bvp4c(odefun,bcfun,solinit,options)
# 说明
sol = bvp4c(odefun,bcfun,solinit) 使用 bcfun 描述的边界条件和初始解估计值 solinit 对 odefun 指定的形式为 y′ = f(x,y) 的微分方程组进行积分。使用 bvpinit 函数创建初始估计值 solinit,该函数还定义了要强制应用 bcfun 中边界条件的点。 示例
sol = bvp4c(odefun,bcfun,solinit,options) 还使用由 options(使用 bvpset 函数创建的参数)定义的积分设置。例如,使用 abstol 和 reltol 选项指定绝对误差容限和相对误差容限,或者使用 fjacobian 选项提供 odefun 的解析偏导数。
# 示例
求解二阶 BVP
使用函数在 Syslab 中求解二阶 BVP。对于此示例,请使用二阶方程
该方程在区间
要在 Syslab 中求解此方程,您需要编写一个将方程表示为一阶方程组的函数、一个边界条件函数和一个初始估计值函数。然后 BVP 求解器使用这三个输入来求解方程。
编写方程代码
编写一个用于编写方程代码的函数。使用代换法
对应的函数是
using TyMath
using TyPlot
function bvpfcn(x, y)
return [y[2]; -y[1]]
end
注意:所有函数都作为局部函数包含在示例末尾。
编写边界条件代码
编写一个函数,以
对应的函数是
function bcfcn(ya, yb)
return [ya[1]; yb[1] - 2]
end
创建初始估计值
使用 bvpinit 函数为方程的解创建初始估计值。由于该方程将
初始估计值的函数接受
function guess(x)
return [sin(x); cos(x)]
end
xmesh = LinRange(0, pi / 2, 5)
solinit = bvpinit(xmesh, guess)
求解方程
结合使用 bvp4c 与导函数、边界条件函数和初始估计值来求解问题。
sol = bvp4c(bvpfcn, bcfcn, solinit)
对解进行绘图
plot(sol.x,sol.y',"-o")
局部函数
此处列出了 bvp4c 用于求解方程的局部函数。
function bvpfcn(x, y) # equation to solve
return [y[2]; -y[1]]
end
# --------------------------------
function bcfcn(ya, yb) # boundary conditions
return [ya[1]; yb[1] - 2]
end
# --------------------------------
function guess(x) # initial guess for y and y'
return [sin(x); cos(x)]
end
# --------------------------------
# 输入参数
odefun — 要求解的函数函数
要求解的函数,指定为定义要积分的函数的函数。odefun 和 bcfun 必须接受相同数量的输入参数。
要编写 odefun 代码,请对标量 x 和列向量 y 使用函数签名 dydx = odefun(x,y)。返回值 dydt 是列向量,该列向量对应于 f(x,y)。odefun 必须同时接受输入参数 x 和 y,即使其中一个参数未在函数中使用也是如此。
例如,要解算
function odefun(t,y)
return 5*y-3;
end
对于方程组,odefun 的输出为向量。向量中的每个元素是一个方程的解。例如,要求解
使用函数:
function odefun(t,y)
return [y[1]+2*y[2]; 3*y[1]+2*y[2]]
end
bvp4c 也可以求解解具有奇异性或具有多点边界条件的问题。
示例: sol = bvp4c(odefun, bcfun, solinit)
未知参数
如果正在求解的 BVP 包含未知参数,您可以使用函数 odefun(x,y,p),其中 p 是参数值向量。当您使用此函数时,BVP 求解器从 solinit 中提供的参数值的初始估计值开始计算未知参数的值。
数据类型: Function
bcfun — 边界条件函数
边界条件,指定为计算边界条件中残差的函数。odefun 和 bcfun 必须接受相同数量的输入参数。
要编写 bcfun 代码,请对列向量 ya 和 yb 使用函数 bcfun(ya,yb)。返回值 res 是列向量,它对应于边界点处边界条件的残差值。
例如,如果 y[a] = 1 且 y[b] = 0,则边界条件函数为
function bcfun(ya,yb)
return [ya[1]-1; yb[1]]
end
由于 y[a] = 1,ya[1]-1 在 x = a 点的残差值应为 0。同样,由于 y[b] = 0,yb[1] 在 x = b 点的残差值应为 0。 强制应用边界条件的边界点 x = a 和 x = b 在初始估计值结构体 solinit 中定义。对于两点边界值问题,a = solinit.x[1] 且 b = solinit.x[end]。
示例: sol = bvp4c(odefun, bcfun, solinit)
未知参数
如果正在求解的 BVP 包含未知参数,您可以使用函数 bcfun(ya,yb,p),其中 p 是参数值向量。当您使用此函数时,BVP 求解器从 solinit 中提供的参数值的初始估计值开始计算未知参数的值。
数据类型: Function
solinit — 解的初始估计值结构体
解的初始估计值,指定为结构体。使用 bvpinit 创建 solinit。
与初始值问题不同,边界值问题可以无解、有有限数量的解或者有无限数量的解。求解 BVP 过程的重要部分是提供所需解的估计值。估计值的准确与否对求解器性能甚至是能否成功计算来说至关重要。有关创建良好初始估计值的一些指导原则,请参阅解的初始估计值。
示例: sol = bvp4c(odefun, bcfun, solinit)
数据类型: Struct
options — options 结构体结构体
options 结构体。使用 bvpset 函数创建或修改 options 结构体。
示例: options = bvpset(;reltol=1e-5,stats=:on) 指定相对误差容限为 1e-5,并打开求解器统计信息的显示。
数据类型: Struct
# 输出参数
sol — 解结构体结构体
解结构体。您可以使用点索引来访问 sol 中的字段,例如 sol.field1。解 (sol.x,sol.y) 在初始网格 solinit.x 中定义的积分区间上是连续的,且在其中有一个连续一阶导数。您可以使用 sol 和 deval 函数来计算该区间中其他点的解。
结构体 sol 包含以下字段。
| 字段 | 描述 |
|---|---|
| x | bvp4c 选择的网格。此网格通常包含不同于初始网格 solinit.x 的点。 |
| y | sol.x 网格点处的 y(x) 近似值。 |
| yp | sol.x 网格点处的 y′(x) 近似值。 |
| parameters | solinit.parameters 中指定的未知参数的最终值。 |
| solver | "bvp4c" |
| stats | 与解相关的计算成本统计信息:网格点的数量、残差以及调用 odefun 和 bcfun 的次数。 |
# 详细信息
多点边界值问题
对于多点边界值问题,会在积分区间内的几个点上强制应用边界条件。
bvp4c 可以求解多点边界值问题(其中 a = a0 < a1 < a2 < ...< an = b 是区间 [a,b] 中的边界点)。点 a1,a2,...,an−1 表示将 [a,b] 分为多个区域的界点。bvp4c 使用从 1 开始的索引按从左到右的顺序(从 a 到 b)枚举区域。在区域 k [ak−1,ak] 中,bvp4c 按如下所示计算导数
在边界条件函数 bcfun(yleft,yright) 中,yleft[:,k] 是 [ak−1,ak] 的左边界的解。同样,yright[:,k) 是区域 k 右边界的解。具体来说是 yleft[:,1] = y(a) 和 yright(:,end] = y(b)。
当您使用 bvpinit 创建初始估计值时,对于每个交界点,请在 xinit 中使用双重项。有关详细信息,请参阅 bvpinit 的参考页。
如果 yinit 为函数,bvpinit 调用 y = yinit(x,k) 以获取区域 k 中的 x 处的初始估计解。在 bpv4c 返回的解结构体 sol 中,sol.x 为每个交界点指定了双重项。sol.y 的对应列分别包含界点的左解和右解。
有关求解三点边界值问题的示例,请参阅求解具有多边界条件的 BVP。
奇异边界值问题
bvp4c 对一类奇异边界值问题求解,其中包括具有未知参数 p 的形式如下的问题
需要区间为 [0, b] 并且 b > 0。当计算偏微分方程 (PDE) 因柱状或球面对称性而产生的 ODE 的平滑解时,通常会出现此类问题。对于奇异问题,应将(常量)矩阵 S 指定为 bvpset 的 singularterm 选项的值,并且 odefun 仅计算 f(x,y,p)。边界条件和初始估计值必须与平滑度 S·y(0) = 0 的必要条件一致。
有关求解奇异边界值问题的示例,请参阅求解具有奇异项的 BVP。
# 算法
bvp4c 是一个有限差分代码,此代码实现 3 阶段 Lobatto IIIa 公式。这是配置公式,并且配置多项式会提供在整个积分区间中处于四阶精度的 C1 连续解。网格选择和误差控制均基于连续解的残差。
配置方法使用点网格将积分区间分为子区间。通过对源于边界条件以及所有子区间上配置条件的线性代数方程全局组求解,求解器会确定数值解。然后,求解器会估计每个子区间上数值解的误差。如果解不满足公差标准,则求解器会调整网格并重复计算过程。您必须提供初始网格以及网格点处解的初始近似估计。