0

我的代码在我的解决方案向量中到处都是零，但我不知道为什么。我已将耦合的二阶 ODE 分解为 4 个一阶 ODE。

我的函数定义为 xp.m

function zprime = f(t,z)
a = 1;
b = 1;
c = 1.5;

zprime = zeros(4,1);

zprime(1) = z(2);
zprime(2) = -a*z(1) + b*(z(3) - z(1));
zprime(3) = z(4);
zprime(4) = -c*(z(3) - z(1));
end

我使用以下命令在 matlab 中运行它：

[t,z] = ode45('xp',[1,100],[0 0 0 0]);

因为我的初始条件都是 0。是我的初始条件给出了 0 解决方案还是其他什么？当我更改 ic 时，解决方案会按预期更改。

谢谢

4

2 回答 2

4

对于您的特定情况，使用 ICs z_0 = [0,0,0,0]，解决方案应该是稳定的，值为z_out = [0,0,0,0]. z = z_0当您插入并通过您的 ODE 求解器运行它时，只需看看您的函数。

zprime(1) = z(2);                       % --->  0
zprime(2) = -a*z(1) + b*(z(3) - z(1));  % ---> -a*(0) + b*(0)
zprime(3) = z(4);                       % ---> 0
zprime(4) = -c*(z(3) - z(1));           % ---> -c*(0-0) = 0

请记住数值解的一般前提。您采用初始条件，并将其输入到导数公式中。这告诉您正在求解的函数的斜率。您使用该斜率来确定函数在未来某个时间（或附近位置或类似位置）的值，然后将其反馈到导数公式中并重新开始。

不同方法（前向/后向 Euler、RK4、...）之间的唯一主要区别是用于确定当前位置斜率的方法。

于 2013-10-17T20:19:13.030 回答
2

您的实际问题与数学有关，而不是 Matlab 或编程。如果你在你的函数中插入零f，你会立即看到它除了零之外没有其他答案。您应该查找平衡点或固定点。即使平衡是不稳定的（想象山顶），恰好在该点没有干扰（外部输入、数值误差）的状态将永远保持平衡。如果您要使用微分方程，最好知道如何找到平衡点，然后如何计算在这些点处评估的雅可比矩阵以确定系统属性。如果您在这方面还有其他问题，我建议您访问math.stackexchange.com。

此外，您正在使用旧方案来调用您的集成函数。你也可以传入你的参数。将您的集成函数定义为主函数中的子函数或单独的函数文件（与函数名相同的文件名——您需要的不是f）

function zprime = f(t,z,a,b,c)
zprime(1,1) = z(2);
zprime(2,1) = -a*z(1) + b*(z(3) - z(1));
zprime(3,1) = z(4);
zprime(4,1) = -c*(z(3) - z(1));

然后，在您的主要功能中，调用

a = 1;
b = 1;
c = 1.5;
[t,z] = ode45(@(t,z)f(t,z,a,b,c),[1 100],[0 0 0 0]);
于 2013-10-17T20:22:55.310 回答