Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question
    
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

1 Answer 1

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.