Matlab：如何在混合初始和终端条件的情况下数值求解ode系统？

Question

我正在尝试使用ode45来解决一组常微分方程组：

[X,Y]=  ode45(@sys,[0, T],y0);

其中，

function dy = sys(t,y)

        dy(1) = f_1(y)
        dy(2) = f_2(y)
        dy(3) = f_3(y)
end

问题在于函数ode45要求y0是初始值[y_1(0), y_2(0), y_3(0)]，但在我的系统中，我只有可用的值[y_2(0), y_3(0), y_3(T)]。

从数学上讲，这组初始/终端条件应该足以确定系统，但我能否通过ode45或MATLAB中的其他函数来处理它呢？

谢谢！

- Vokram

---

我对这个问题非常感兴趣，但恐怕无法提供帮助；我以前从未遇到过这种问题... 我知道 ode45 可以向后积分（只需使用 tspan = [tend tstart]），因此您可以构造一个迭代方案来获得 y_1，使得满足 y_3(0) 和 y_3(T)。不用说，这可能会非常缓慢和笨拙，但它是一种解决方案。我会关注这个问题 :) 您能发布方程式和初始/终端条件吗？ - Rody Oldenhuis

@RodyOldenhuis 我想我找到了解决这个问题的方法。 - Vokram

Answer 1

在研究了一下Matlab文档后，我认为更优雅的方法是使用bvp4c函数。bvp4c是一个专门设计用于处理边界值问题的函数，而不是像ode**那样仅适用于初值问题。实际上，在Matlab中还有一整套其他函数，如deval和bvpinit，这些函数真正方便了bvp4c的使用。这是Matlab文档的链接：http://www.mathworks.com/help/matlab/math/boundary-value-problems.html。

这里我将发布一个简短（也许有点牵强）的示例：

function [y1, y2, y3] = test(start,T)

solinit = bvpinit(linspace(0,3,10), [1,1,0]);
sol = bvp4c(@odefun,@bvpbc,solinit);

tspan = linspace(start,T,100);
S = deval(sol, tspan);
y1 = S(1,:);
y2 = S(2,:);
y3 = S(3,:);

plot (tspan,y1)

figure
plot (tspan,y2)

figure
plot (tspan,y3)


%% system definition & BVCs

function dydx = odefun(t,y)
    dydx(1) = y(1) + y(2) + t;

    dydx(2) = 2*y(1) + y(2);

    dydx(3) = 3 * y(1) - y(2);
end

function res = bvpbc(y0,yT)
   res= [y0(3) yT(2) yT(3)];
end

end

test函数会分别输出3个解向量y1、y2和y3，并绘制它们的图形。下面是Matlab绘制的变量路径： 此外，我发现WMU的Jake Blanchard教授的这个视频非常有帮助。

- Vokram

---

+2：看起来你解决了它！我今天学到了新东西，谢谢 :) - Rody Oldenhuis

+2 分钟，因为 @RodyOldenhuis 不知道 BVP 求解器而感到惊讶！ :-) 如果您喜欢，可以接受自己的答案。 - horchler

Answer 2

一种方法是使用射击法迭代求解未知的初始状态y_1(0)，以满足所需的最终状态y_3(T)。

迭代过程通过求解非线性方程F = 0进行：

F(y_1(0)) = Y_3(T) - y_3(T)

其中，大写函数Y是通过从一组初始条件向时间进行正向ODE积分获得的解。任务是猜测y_1(0)的值以获得F=0。

由于我们现在正在解决一个非线性方程，因此所有通常的方法都适用。具体来说，您可以使用二分法或牛顿法更新未知初始状态y_1(0)的猜测。请注意以下几点：

1. ODE在每次非线性求解迭代中在[0,T]上进行积分。
2. 对于F=0可能会有多个解，这取决于您的ODE结构。
3. 牛顿法可能比二分法收敛更快，但如果您不能提供y_1(0)的良好起始猜测，则可能在数值上不稳定。

使用现有的MATLAB函数，有界标量非线性求解器FMINBND可能是作为非线性求解器的一个不错选择。

希望这可以帮助到您。

- Darren Engwirda

---

+1：这本来是我会采取的方法，但我认为 @Vokram 的答案应该被投票置顶 :) - Rody Oldenhuis