matlab - 二阶系统的中点欧拉方法，不能让它工作

Question

我需要绘制这个简单的系统：

x'' = -x

使用中点欧拉。

u1 = -x , u2 = -x'

u1' = u2
u2' = -x

u1(n+1) = u1(n) + h*？
u2(n+1) = u2(n) + h*f((1/2)*(u1(n) + u1(n+1))

我们不知道 u1(n+1)。我尝试用 u1(n+1) = u1(n) + h*u2(n) 来近似它

u2(1+i) = u2(i) + h * ((-1/2) * (u1(n) + u1(n+1))

然后我们有 u2(i+1) 和 u2(i)。中点值为 (u2(1+i) - u(i))/2

u1(q+i) = u1(i) + h*中点

当我绘制它时，结果是一些可怕的发散线，而不是振荡函数。怎么了？

clear all, close all, clc

h = 0.0005; % steplength
nos = 500000; % number of steps

x = zeros(1,nos);
xp = zeros(1,nos);
energy=zeros(1,nos);

% Starting positions
x(1) = 1;
xp(1) = 0;

for i=1:nos-1    



    xpp = -x(i);

    xTAYLOR = x(i) + h*xp(i);

    xp(1+i) = x(i) + (-1*((1/2)*(x(i) + xTAYLOR)));

    xpHALF = (xp(1+i) - x(i))/2;

    x(1+i) = x(i) + h*xpHALF;


end

plot(x)

Answer 1

好的，也许我被击中了，但我看到问题的提出方式存在一个小问题。

u1 = -x , u2 = -x'

u1' = u2
u2' = -x

让我试试这个：

u1 = -x
u2 = -x'

u1' = -x' = u2
u2' = -x'' = x = -u1

当我这样做时，我会得到一些甜蜜的摆动！这是我的尝试，只是使用ode45，但是您应该对符号更改有类似的运气。

[t2,y2] = ode45(@(t,x) [x(2) -x(1)]',[0 250],[0 1]);

嗯！！