计算科学 - 在 MatLab 中求解 BVP 时如何实现积分条件 - 吾爱随笔录

在 MatLab 中求解 BVP 时如何实现积分条件

计算科学 matlab 颂边界条件

2021-12-26 17:37:57

我正在尝试使用 MATLAB 的bvp4c函数来解决 ODE 的耦合系统。我想强加其中是第一个状态变量。

\int_{0}^{π} y_{1} (t) y_{1} (t) d t = 1,

$\int_{0}^{\pi} y_{1}(t) y_{1}(t) dt = 1,$

y_{1}

$y_{1}$

另一篇文章建议在 ODE 系统中添加一个额外的方程。即

I (t) = \int_{0}^{t} y_{1} (s) y_{1} (s) d s .

$I(t) = \int_{0}^{t} y_{1}(s) y_{1}(s) ds \, .$

然后附加的边界条件将强制执行我想要的积分条件。但是，我不确定如何在 MATLAB 中编写代码。谁能帮我这个？ $I(\pi)-1=0$

按照尼科瓜罗的评论，积分方程可以写为

\frac{d}{d t} I (t) = y_{1} (t) y_{1} (t) .

$\frac{d}{dt}I(t) = y_{1}(t) y_{1}(t) \, .$

相应地更新代码，我有

这是我正在使用的整个代码bvp4c（如基里尔说，包含整个代码会很有用）：

s = 100;
solinit = bvpinit(linspace(0,pi,1001),@mat4init,s);
sol = bvp4c(@odeproblem,@eightbc,solinit);

x = linspace(0,pi,1001);
y = deval(sol,x);
trapz(y(1,:).*y(1,:))
y_bar = y(1,:)./norm(y(1,:),2);
plot(x,y(1,:))

fprintf('The intial guess for s was %d. The value solved for was %d.\n',s,sol.parameters)

function dydt = odeproblem(t,y,s)
b = pi;
dydt(1) = y(2);
dydt(2) = y(3);
dydt(3) = y(4);
dydt(4) = y(5);
dydt(5) = y(6);
dydt(6) = y(7);
dydt(7) = y(8);
dydt(8) = (-b^4*y(1) + 2*b^2*y(3) - (b^4-2)*y(5) - 2*s^2*b^2*y(7))/s^2;
dydt(9) = y(1)*y(1);
end

function res = eightbc(ya,yb,s)
res = [ya(1); yb(1); ya(3); yb(3); ya(5); yb(5); ya(7); yb(7); ya(9); yb(9)-1];
end

function yinit = mat4init(x)
yinit = [  cos(2*x), -2*sin(2*x), -4*cos(2*x), 8*sin(2*x),... 
    16*cos(2*x),  -32*sin(2*x), -64*cos(2*x), 128*sin(2*x), 1];
end

这是解决方案的图：

所有其他边界条件似乎都得到满足（我没有包括显示的图）。但积分条件仍然不是。

1个回答

注释详细说明了通过将积分定义为 ODE 的另一部分并在其上添加边界条件来解决此问题的一种完全有效的方法。我将讨论解决此问题的另一种方法。

bvp4c是一个两点边值问题求解器。您需要能够指定使用多个点的条件，或者更好的是，使用当前解决方案的连续扩展，而不是两点边界值问题求解器。您可以使用DifferentialEquations.jl 执行类似的操作，因为它允许您定义具有解的连续扩展的残差。

function bc3(residual, sol)
    residual[1] = quadgk((t)->sol(t;idxs=1)^2,0,pi)
    residual[2:8] .= 0 #there's no other boundary conditions mentioned?
end

这是有效的，因为它是使用顺序匹配插值sol(t)的解决方案。t

好的，所以可以直接处理这个，但有什么区别？这是更有趣的部分。两种不同的情况。当使用射击型边值问题求解器时，自然而然地进行连续扩展，并且对初始条件进行求根，将 BVP 视为具有未知初始条件的 IVP。因此，使用这种连续扩展条件没有缺点（只要 BVP 求解器允许）。但拍摄方法当然有一个缺点，即它们不能用于对初始条件敏感或以奇点开始的问题。

这不适用于搭配或 MIRK（隐式 RK）类型的方法。在这些情况下，一个网格空间并通过 rootfinder中的每个点放宽解。在这种情况下，“雅可比”是每个点对每个其他点的依赖关系。由于对于 Runge-Kutta 或搭配方法，您只使用之前的值来计算新值，因此这通常具有带状结构。但是，这是假设边界仅取决于和。如果你施加一个依赖于空间中所有点的边界条件，比如积分的这种正交形式，那么整个系统的雅可比行列就有一个完整的列，它破坏了带状结构。 $x_i$ $y(x_i)$ $y(x_1)$ $y(x_{end})$

由于雅可比矩阵的大小M*N*S是M空间点N的数量，是系统中 ODE 的数量，并且S是 ODE 中的阶段数，即使对于简单的问题，这也可能非常大。但如前所述，两点 BVP 问题使依赖关系成为S带状矩阵，因此很容易解决。bvp4c是许多此类求解器之一，因此采用两点 BVP 形式，因为它可以大大简化求解该系统。

如果不是两点 BVP，则此 Jacobian 矩阵是稀疏矩阵。要么它实际上存储为稀疏矩阵，要么如果问题足够小，则使用密集矩阵，但反演算法不能使用带状形式，因此采用更慢的\. 因此，只有在必要时才应使用像这样泛化的求解器。

让我总结一下：

如果您可以进行转换以不具有更复杂的边界条件，那么这可能是您最好的选择，如果它不难做到的话。积分可以作为额外的部分添加到 ODE 系统中，因为 ODE 求解器实际上只是一种求积方法。
如果问题比较困难并且不可能进行这样的转换，那么射击型方法是一个不错的选择，并且不会因遇到这些情况而产生开销（您只需要找到一个 BVP 求解器即可它）。
射击方法要求问题对初始条件不敏感或具有初始奇点。如果不是这种情况，则需要使用搭配或 MIRK 类型的方法。然而，这些方法中的大多数，例如bvp4c，不允许这种形式，因为它们可以通过假设问题具有两点条件来极大地优化。少数可以处理这种情况的求解器只能作为最后的手段。

希望这些解释有助于让您清楚，评论中的解决方法不仅仅是完成它的技巧，而且可能是更好的方法。

其它你可能感兴趣的问题

上一篇使用 Crank-Nicolson 方法求解空间和时间方程下一篇用于并行化有限元模拟的简单快速图形聚类