bvp5c
求解边界值问题 - 五阶方法
说明
示例
求解二阶 BVP
使用函数在 MATLAB® 中求解二阶 BVP。对于此示例,请使用二阶方程
.
该方程在区间 上定义并受限于边界条件
,
.
要在 MATLAB 中求解此方程,您需要编写一个将方程表示为一阶方程组的函数、一个边界条件函数和一个初始估计值函数。然后 BVP 求解器使用这三个输入来求解方程。
编写方程代码
编写一个用于编写方程代码的函数。使用代换法 和 将方程重写为一阶方程组。
,
.
对应的函数是
function dydx = bvpfcn(x,y) dydx = zeros(2,1); dydx = [y(2) -y(1)]; end
注意:所有函数都作为局部函数包含在示例末尾。
编写边界条件代码
编写一个函数,以 形式编写边界条件代码。在此形式中,边界条件是
,
.
对应的函数是
function res = bcfcn(ya,yb) res = [ya(1) yb(1)-2]; end
创建初始估计值
使用 bvpinit
函数为方程的解创建初始估计值。由于该方程将 与 联系起来,合理的估计是该解涉及三角函数。使用由积分区间的五个点构成的网格。网格中的第一个和最后一个值是求解器应用边界条件的位置。
初始估计值的函数接受 作为输入,并返回对 和 值的估计值。函数是
function g = guess(x) g = [sin(x) cos(x)]; end
xmesh = linspace(0,pi/2,5); solinit = bvpinit(xmesh, @guess);
求解方程
结合使用 bvp5c
与导函数、边界条件函数和初始估计值来求解问题。
sol = bvp5c(@bvpfcn, @bcfcn, solinit);
对解进行绘图
plot(sol.x, sol.y, '-o')
局部函数
此处列出了 bvp5c
用于求解方程的局部函数。
function dydx = bvpfcn(x,y) % equation to solve dydx = zeros(2,1); dydx = [y(2) -y(1)]; end %-------------------------------- function res = bcfcn(ya,yb) % boundary conditions res = [ya(1) yb(1)-2]; end %-------------------------------- function g = guess(x) % initial guess for y and y' g = [sin(x) cos(x)]; end %--------------------------------
比较 bvp4c
和 bvp5c
求解器
用两个不同求解器以宽松误差容限求解 BVP,并比较结果。
假设使用二阶 ODE
.
该方程在区间 上定义并受限于边界条件
,
.
要在 MATLAB® 中求解此方程,您需要编写一个将方程表示为一阶方程组的函数,编写一个边界条件函数,设置一些选项值,并创建初始估计值。然后 BVP 求解器使用这四个输入来求解方程。
编写方程代码
使用代换法 和 ,您可以将 ODE 重写为一阶方程组
,
.
对应的函数是
function dydx = bvpfcn(x,y) dydx = [y(2) -2*y(2)/x - y(1)/x^4]; end
注意:所有函数都作为局部函数包含在示例末尾。
编写边界条件代码
边界条件函数要求边界条件的形式为 。在此形式中,边界条件是
,
.
对应的函数是
function res = bcfcn(ya,yb) res = [ya(1) yb(1)-sin(1)]; end
设置选项
使用 bvpset
打开求解器统计信息的显示,并指定宽松误差容限以突出求解器之间误差控制的差异。此外,为了提高效率,请指定解析雅可比矩阵
.
返回雅可比矩阵值的对应函数为
function dfdy = jac(x,y) dfdy = [0 1 -1/x^4 -2/x]; end
opts = bvpset('FJacobian',@jac,'RelTol',0.1,'AbsTol',0.1,'Stats','on');
创建初始估计值
使用 bvpinit
创建解的初始估计值。指定一个常量函数作为初始估计值,初始网格包含区间 中的 10 个点。
xmesh = linspace(1/(3*pi), 1, 10); solinit = bvpinit(xmesh, [1; 1]);
求解方程
用 bvp4c
和 bvp5c
分别求解方程。
sol4c = bvp4c(@bvpfcn, @bcfcn, solinit, opts);
The solution was obtained on a mesh of 9 points. The maximum residual is 9.794e-02. There were 157 calls to the ODE function. There were 28 calls to the BC function.
sol5c = bvp5c(@bvpfcn, @bcfcn, solinit, opts);
The solution was obtained on a mesh of 11 points. The maximum error is 6.742e-02. There were 244 calls to the ODE function. There were 29 calls to the BC function.
绘制结果
绘制 的两次计算结果和解析解以进行比较。解析解是
,
.
xplot = linspace(1/(3*pi),1,200); yplot = [sin(1./xplot); -cos(1./xplot)./xplot.^2]; plot(xplot,yplot(1,:),'k',sol4c.x,sol4c.y(1,:),'r*',sol5c.x,sol5c.y(1,:),'bo') title('Comparison of BVP Solvers with Crude Error Tolerance') legend('True','BVP4C','BVP5C') xlabel('x') ylabel('solution y')
绘图验证了 bvp5c
直接控制计算中的真误差,而 bvp4c
仅间接控制计算中的真误差。在更严格的误差容限下,求解器之间的这种差异没有这么明显。
局部函数
此处列出了 BVP 求解器用于求解问题的局部函数。
function dydx = bvpfcn(x,y) % equation to solve dydx = [y(2) -2*y(2)/x - y(1)/x^4]; end %--------------------------------- function res = bcfcn(ya,yb) % boundary conditions res = [ya(1) yb(1)-sin(1)]; end %--------------------------------- function dfdy = jac(x,y) % analytical jacobian for f dfdy = [0 1 -1/x^4 -2/x]; end %---------------------------------
输入参数
odefun
— 要求解的函数
函数句柄
要求解的函数,指定为定义要积分的函数的函数句柄。odefun
和 bcfun
必须接受相同数量的输入参量。
要编写 odefun
代码,请对标量 x
和列向量 y
使用函数签名 dydx = odefun(x,y)
。返回值 dydx
是数据类型为 single
或 double
的列向量,该列向量对应于 f(x,y)。odefun
必须同时接受输入参量 x
和 y
,即使其中一个参量未在函数中使用也是如此。
例如,要解算 ,请使用此函数:
function dydx = odefun(x,y) dydx = 5*y-3; end
对于方程组,odefun
的输出为向量。向量中的每个元素是一个方程的解。例如,要求解
使用函数:
function dydx = odefun(x,y) dydx = zeros(2,1); dydx(1) = y(1)+2*y(2); dydx(2) = 3*y(1)+2*y(2); end
bvp5c
也可以求解解具有奇异性或具有多点边界条件的问题。
示例: sol = bvp5c(@odefun, @bcfun, solinit)
未知参数
如果正在求解的 BVP 包含未知参数,您可以使用函数签名 dydx = odefun(x,y,p)
,其中 p
是参数值向量。当您使用此函数签名时,BVP 求解器从 solinit
中提供的参数值的初始估计值开始计算未知参数的值。
数据类型: function_handle
bcfun
— 边界条件
函数句柄
边界条件,指定为计算边界条件中残差的函数句柄。odefun
和 bcfun
必须接受相同数量的输入参量。
要编写 bcfun
代码,请对列向量 ya
和 yb
使用函数签名 res = bcfun(ya,yb)
。返回值 res
是数据类型为 single
或 double
的列向量,它对应于边界点处边界条件的残差值。
例如,如果 y(a) = 1 且 y(b) = 0,则边界条件函数为
function res = bcfun(ya,yb) res = [ya(1)-1 yb(1)]; end
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 = bvp5c(@odefun, @bcfun, solinit)
未知参数
如果正在求解的 BVP 包含未知参数,您可以使用函数签名 res = bcfun(ya,yb,p)
,其中 p
是参数值向量。当您使用此函数签名时,BVP 求解器从 solinit
中提供的参数值的初始估计值开始计算未知参数的值。
数据类型: function_handle
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
包含以下字段。
字段 | 描述 |
---|---|
|
|
|
|
|
|
|
|
|
|
| 与解相关的计算成本统计信息:网格点的数量、残差以及调用 |
详细信息
多点边界值问题
对于多点边界值问题,会在积分区间内的几个点上强制应用边界条件。
bvp5c
可以求解多点边界值问题(其中 a = a0 < a1 < a2 < ...< an = b 是区间 [a,b] 中的边界点)。点 a1,a2,...,an−1 表示将 [a,b] 分为多个区域的界点。bvp5c
使用从 1 开始的索引按从左到右的顺序(从 a 到 b)枚举区域。在区域 k [ak−1,ak] 中,bvp5c
按如下所示计算导数
yp = odefun(x,y,k)
在边界条件函数 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。
奇异边界值问题
bvp5c
对一类奇异边界值问题求解,其中包括具有未知参数 p
的形式如下的问题
需要区间为 [0, b] 并且 b > 0。当计算偏微分方程 (PDE) 因柱状或球面对称性而产生的 ODE 的平滑解时,通常会出现此类问题。对于奇异问题,应将(常量)矩阵 S
指定为 bvpset
的 'SingularTerm'
选项的值,并且 odefun
仅计算 f(x,y,p)。边界条件和初始估计值必须与平滑度 S·y(0) = 0 的必要条件一致。
有关求解奇异边界值问题的示例,请参阅求解具有奇异项的 BVP。
算法
bvp5c
是一个有限差分代码,此代码实现 4 阶段 Lobatto IIIa 公式 [1]。这是配置公式,并且配置多项式会提供在整个 [a,b]
中具有一致五阶精度的 CC1 连续解。该公式作为隐式龙格-库塔公式实现。bvp5c
和 bvp4c
之间的一些差异是:
bvp5c
直接求解代数方程。bvp4c
使用解析凝聚。bvp4c
直接处理未知参数。bvp5c
使用未知参数的平凡微分方程扩充方程组。
参考
[1] Shampine, L.F., and J. Kierzenka. "A BVP Solver that Controls Residual and Error." J. Numer. Anal. Ind. Appl. Math. Vol. 3(1-2), 2008, pp. 27–41.
扩展功能
基于线程的环境
使用 MATLAB® backgroundPool
在后台运行代码或使用 Parallel Computing Toolbox™ ThreadPool
加快代码运行速度。
此函数完全支持基于线程的环境。有关详细信息,请参阅在基于线程的环境中运行 MATLAB 函数。
版本历史记录
在 R2006b 中推出
MATLAB 命令
您点击的链接对应于以下 MATLAB 命令:
请在 MATLAB 命令行窗口中直接输入以执行命令。Web 浏览器不支持 MATLAB 命令。
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)