It's a 2D ODE in my assignment,
$\begin{bmatrix}x_1'(t)\\x_2'(t)\end{bmatrix}=\begin{bmatrix}-1&1\\1&-1000\end{bmatrix}\begin{bmatrix}x_1(t)\\x_2(t)\end{bmatrix}+\begin{bmatrix}2\sin(t)\\1000(\cos(t)-\sin(t))\end{bmatrix}$ on $t\in[0,30]$
with initial values
$x_1(0)=1,x_2(0) = 2$
I tried to get the analytical solution in Mathematica, it gives me a long long expression with a sharp & nearly 90 degree descending at about $t_p=0.72$.
Apprently, this is an ill equation. I tried some 3 order LMM (Linear Multi-Step Method) recursions, say
$y_{n+1}=y_{n-2}+\frac{h}{4}\left[3f(x_{n+1},y_{n+1})+9f(x_{n-1},y_{n-1})\right]$
or
$y_{n+1}=y_{n-1}+\frac{h}{3}\left[7f(x_n,y_n)-2f(x_{n-1},y_{n-1})+f(x_{n-2},y_{n-2})\right]$
but they all explodes after $t_p$ (with slightly larger stepsize, even explodes from the beginning!), no matter how small stepsize I take.
I'm really confused what to do next. Could someone enlighten me?
The assignment asks me to use 3-steps 3 order LMM, but besides this, I'm also curious for an ill ODE, is there any general/common tools to deal with this kind of equations?
The first method is stable where $(1-\frac34z)q^3-\frac94zq-1=0$, $z=\lambda h$, has all roots $q$ inside or on the unit circle if $Im(z)\le 0$. For $z=0$ these are the third unit roots, so you get zero-stability. In first order in $z$, the equation is the same as $q^3-\frac94zq-(1+\frac34z)=O(z^2)$. As the constant term is outside the unit circle, at least one root $q$ has to be outside the unit circle, making the method unstable or weakly stable. This means errors in the initialization/bootstrap or floating point errors get magnified exponentially.
For the second method, the same procedure gives the polynomial $z^3-\frac73zq^2-(1-\frac23z)q-\frac13z=0$ with roots $-\frac13z - \frac29z^2-…$, $-1 + \frac53z - \frac5{18}z^2+…$ and $1 + z + \frac12z^2 +…$ The middle one has a real part smaller $-1$ for $Im(z)<0$, so again the method is inherently unstable.
In summary, it is not surprising that you get unstable results independent of the bootstrap implementation, as the error term reaches the level of the solution values at about $t=0.02$ and reaches overflow sometime after $t=0.4$ for the second method.
For the system itself, if one defines the residuals to a guess of the asymptotic functions $u=x-\sin t$, $v=y-\cos t+\sin t$, $(x,y)=(x_1,x_2)$, then $$ u'+u-v =0\\ v'-u+1000v =2\sin t+\cos t $$ So $u$ falls slowly to $v$ with rate $1$, while $v$ falls quickly with rate 1000 to $10^{-3}(u+2\sin t+\cos t)$. Inserting $u\sim v$ gives $v\sim 999·10^{-6}·(2\sin t+\cos t)$. This is, of course, circular and would require a proper eigen-decomposition for exact results, but for an approximate picture of what happens this is sufficient, the residual level is below the resolution of "not zoomed" plots.
Plots of the solution (with the asymptotes) and step sizes