# How do I use MATLAB to solve an Euler Backward problem for a Differential Equation?

Is there a way to use MATLAB to solve an Euler Backward problem when the function I have is a differential equation?

My original problem was to simulate an electromagnetic and an electrical field for a certain interval in R^2 and in this interval an electron would enter with an initial speed in x, y direction. I then plotted the path the electron would take inside the field, but now I'm trying to approximate the path using Euler Backward.

Our differential equation in MATLAB looks like this:

dy(1)=y(3);
dy(2)=y(4);
dy(3)=y(4)+1;
dy(4)=-y(3)


Where y=[x ; y ; u ; v] where u is the derivative of x and v is that of y.

-
Are you required to write your own code to implement backward Euler, or can you use one of Matlab's prewritten routines? –  Chris Taylor May 13 '11 at 11:05
EDIT: Yes, we are required to write our own code to implement backward Euler, we can not use one of Matlab's prewritten routines. –  Adam May 13 '11 at 13:05
Does anybody have any suggestions? –  Adam May 13 '11 at 19:03

Given the vector-valued differential equation $\mathbf{y}'(t) = \mathbf{f}(\mathbf{y}(t))$, the Euler backward scheme is given by discretizing the time derivative and evaluating $\mathbf{f}$ at the value of $\mathbf{y}$ to be determined:

$$\frac{\mathbf{y}_{n+1} - \mathbf{y}_n}{h} = \mathbf{f}(\mathbf{y}_{n+1})$$

where $h$ is the time step and $t = nh$, so you have

$$\mathbf{y}_{n+1} - \mathbf{y}_n - h\mathbf{f}(\mathbf{y}_{n+1}) = 0$$

which is an implicit equation for $\mathbf{y}_{n+1}$. There are two interesting facts about the case you're trying to solve:

(1) the equation is basically linear, i.e. $\mathbf{f}(\mathbf{y}) = \mathbf{b} + A\mathbf{y}$ where $\mathbf{b}=(0,0,1,0)$ and $A$ is a matrix.

(2) the differential equation for $y_3$ and $y_4$ has no dependence on $y_1$ and $y_2$, implying that you can solve for $y_3$ and $y_4$ separately, and then calculate $y_1$ and $y_2$.

If we consider that $\mathbf{f}$ is linear, then we get the equation

$$\mathbf{y}_{n+1} - \mathbf{y}_n - hA\mathbf{y}_{n+1} - h\mathbf{b} = 0$$

which can be solved explicitly for $\mathbf{y}_{n+1}$, giving

$$\mathbf{y}_{n+1} = (I - hA)^{-1}(\mathbf{y}_n + h\mathbf{b})$$

where $I$ is the identity matrix. This requires finding the inverse of a 4 x 4 matrix which, whilst it is entirely possible using a computer, isn't necessarily desirable for stability reasons. It would be better to just consider the last two components of the equation, the $y_3$ and $y_4$ components. After manually inverting the 2 x 2 matrix you get

$$y_3^{n+1} = \frac{y_3^n + h y_4^n + h}{1+h^2}$$

$$y_4^{n+1} = \frac{y_4^n - hy_3^n}{1+h^2}$$

which you should easily be able to program in Matlab, and then use the values for $y_3$ and $y_4$ to calculate $y_2$ and $y_1$.

Note that I've given you a very specific implementation of backward Euler which is tailored for the problem at hand. A more general backward Euler solver would have to use a root-finding technique such as Newton-Raphson iteration to calculate $\mathbf{y}_{n+1}$ at each step.

-