I want to consider this differential system: $$ \ \frac{dx}{dt} = -y(t)\\ \frac{dy}{dt} = \ x(t) $$ where $t>0$ with initial condition$ (x(0),y(0))=(1,0).$
First I want to show that this differential equation admits an invariant of $I = x(t)^2 + y(t)^2$ Also, Can someone help me to figure out if Forward Euler, backwards Euler, or implicit trapezoidal rule admit invariants similar to $I = x(t)^2 + y(t)^2$?
I wonder how I found this old question. The right way to consider this is interpreting your equation (which is the harmonic oscillator) as a Hamiltonian system with Hamiltonian $H= (x^2+y^2)/2$.
Then, the way to find an invariant is to look for a so-called modified Hamiltonian which is solved exactly by the method and hence an invariant.
Good candidates for such methods go by the name symplectic integrators since their modified equations are Hamiltonian. None of your methods is symplectic and you will probably have a hard time finding invariants in the modified equations (which are then not Hamiltonian) ;) Try implicit midpoint & symplectic Euler, the latter has the invariant $H=x^2+y^2+hxy$ where $h$ is your time-step.
I would like to explain this, using explicit Euler, because it is the most elementary. We have: \begin{align} \begin{pmatrix} x \\ y \end{pmatrix}' = \begin{pmatrix} -y \\ x \end{pmatrix} = f(x,y) \end{align} Explicit euler looks like this: $z_{k+1} = z_k +h \cdot f(z_k)$, applied to our problem it becomes: \begin{align} \begin{pmatrix} x_{k+1} \\ y_{k+1} \end{pmatrix} = \begin{pmatrix} x_{k} \\ y_{k} \end{pmatrix} +h \begin{pmatrix} -y_{k} \\ x_{k} \end{pmatrix} = \begin{pmatrix} x_{k}-hy_k \\ y_{k}+h x_k \end{pmatrix} \end{align} So if we want to know, whether $x(t)^2+y(t)^2$ are constant, we check \begin{align} x_{k+1}^2+y_{k+1}^2 = (x_k -h y_k)^2 +(y_k+ h x_k)^2 = (x_k^2+y_k^2)(1+h^2) \end{align} Give the initial values $x_0=1$ and $y_0=0$ and using recursion, we get \begin{align} x_{k+1}^2+y_{k+1}^2 &=(x_k^2+y_k^2)(1+h^2) =(x_{k-1}^2+y_{k-1}^2)(1+h^2)^2 =\ ... \ =(x_0^2+y_0^2)(1+h^2)^{k+1}\\ &= (1+h^2)^{k+1} \end{align} So it is no invariant.
Next implicit euler. The method reads $z_{k+1} = z_k +h \cdot f(z_{k+1})$ So we get \begin{align} & \begin{pmatrix} x_{k+1} \\ y_{k+1} \end{pmatrix} = \begin{pmatrix} x_{k} \\ y_{k} \end{pmatrix} +h \begin{pmatrix} -y_{k+1} \\ x_{k+1} \end{pmatrix} = \begin{pmatrix} x_{k}-hy_{k+1} \\ y_{k}+h x_{k+1} \end{pmatrix} \\ & \Leftrightarrow \begin{pmatrix} x_{k+1}+ h y_{k+1} \\ y_{k+1} -h x_{k+1} \end{pmatrix} = \begin{pmatrix} x_{k} \\ y_{k} \end{pmatrix} \\ & \Leftrightarrow \begin{pmatrix} 1 & h \\ 1 & -h \end{pmatrix} \begin{pmatrix} x_{k+1} \\ y_{k+1} \end{pmatrix} = \begin{pmatrix} x_{k} \\ y_{k} \end{pmatrix} \end{align} So as you can see, you end up with a linear system you have to solve for $x_{k+1} $ and $y_{k+1}$. This is the price you pay for an implicit scheme. If you solve for $x_{k+1} $ and $y_{k+1}$ you can again check the invariant. Good luck,
Thomas
To a given problem $ z' = f(z) $, given an inital value $z(0)=z_0$, we can approximate the solution of this ODE (ordinary differential equation), by evaluating the our approximation at discrete points. Thus, the FE can be written as $z_{k+1} = z_k +h \cdot f(z_k)$ with som stepsize h. The smaller h, the more detailed you approximate the solution. For backward Euler the methods reads $z_{k+1} = z_k +h \cdot f(z_{k+1})$.
– Thomas Nov 12 '12 at 21:44
