MATLAB code help. Backward Euler method

Question

Here is the MATLAB/FreeMat code I got to solve an ODE numerically using the backward Euler method. However, the results are inconsistent with my textbook results, and sometimes even ridiculously inconsistent. What is wrong with the code?

function [x,y] = backEuler(f,xinit,yinit,xfinal,h)

    %f - this is your y prime
    %xinit - initial X
    %yinit - initial Y
    %xfinal - final X
    %h - step size

    n = (xfinal-xinit)/h; %Calculate steps

    %Inititialize arrays...
    %The first elements take xinit and yinit corespondigly, the rest fill with 0s.
    x = [xinit zeros(1,n)];
    y = [yinit zeros(1,n)];

    %Numeric routine
    for i = 1:n
        x(i+1) = x(i)+h;
        ynew = y(i)+h*(f(x(i),y(i)));
        y(i+1) = y(i)+h*f(x(i+1),ynew);
    end
end

Answer 1

Your method is a method of a new kind. It is neither backward nor forward Euler. :-)

Forward Euler: y1 = y0 + h*f(x0,y0)

Backward Euler solve in y1: y1 - h*f(x1,y1) = y0

Your method: y1 = y0 +h*f(x0,x0+h*f(x0,y0))

Your method is not backward Euler.

You don't solve in y1, you just estimate y1 with the forward Euler method. I don't want to pursue the analysis of your method, but I believe it will behave poorly indeed, even compared with forward Euler, since you evaluate the function f at the wrong point.

Here is the closest method to your method that I can think of, explicit as well, which should give much better results. It's Heun's Method:

y1 = y0 + h/2*(f(x0,y0) + f(x1,x0+h*f(x0,y0)))

Answer 2

The only issue I can spot is that the line:

n=(xfinal-xinit)/h

Should be:

n = abs((xfinal-xinit)/h)

To avoid a negative step count. If you are moving in the negative-x direction, make sure to give the function a negative step size.

Your answers probably deviate because of how coarsely you are approximating your answer. To get a semi-accurate result, deltaX has to be very very small and your step size has to be very very very small.

PS. This isn't the "backward Euler method," it is just regular old Euler's method.

If this is homework please tag it so.

Answer 3

Have a look at numerical recipes, specifically chapter 16, integration of ordinary differential equations. Euler's method is known to have problems:

There are several reasons that Euler’s method is not recommended for practical use, among them, (i) the method is not very accurate when compared to other, fancier, methods run at the equivalent stepsize, and (ii) neither is it very stable

So unless you know your textbook is using Euler's method, I wouldn't expect the results to match. Even if it is, you probably have to use an identical step size to get an identical result.

Answer 4

Unless you really want to solve an ODE via Euler's method that you've written by yourself you should have a look at built-in ODE solvers.

On a sidenote: you don't need to create x(i) inside the loop like this: x(i+1) = x(i)+h;. Instead, you can simply write x = xinit:h:xfinal;. Also, you may want to write n = round(xfinal-xinit)/h); to avoid warnings.

---

Here are the solvers implemented by MATLAB.

ode45 is based on an explicit Runge-Kutta (4,5) formula, the Dormand-Prince pair. It is a one-step solver – in computing y(tn), it needs only the solution at the immediately preceding time point, y(tn-1). In general, ode45 is the best function to apply as a first try for most problems.

ode23 is an implementation of an explicit Runge-Kutta (2,3) pair of Bogacki and Shampine. It may be more efficient than ode45 at crude tolerances and in the presence of moderate stiffness. Like ode45, ode23 is a one-step solver.

ode113 is a variable order Adams-Bashforth-Moulton PECE solver. It may be more efficient than ode45 at stringent tolerances and when the ODE file function is particularly expensive to evaluate. ode113 is a multistep solver — it normally needs the solutions at several preceding time points to compute the current solution.

The above algorithms are intended to solve nonstiff systems. If they appear to be unduly slow, try using one of the stiff solvers below.

ode15s is a variable order solver based on the numerical differentiation formulas (NDFs). Optionally, it uses the backward differentiation formulas (BDFs, also known as Gear's method) that are usually less efficient. Like ode113, ode15s is a multistep solver. Try ode15s when ode45 fails, or is very inefficient, and you suspect that the problem is stiff, or when solving a differential-algebraic problem.

ode23s is based on a modified Rosenbrock formula of order 2. Because it is a one-step solver, it may be more efficient than ode15s at crude tolerances. It can solve some kinds of stiff problems for which ode15s is not effective.

ode23t is an implementation of the trapezoidal rule using a "free" interpolant. Use this solver if the problem is only moderately stiff and you need a solution without numerical damping. ode23t can solve DAEs.

ode23tb is an implementation of TR-BDF2, an implicit Runge-Kutta formula with a first stage that is a trapezoidal rule step and a second stage that is a backward differentiation formula of order two. By construction, the same iteration matrix is used in evaluating both stages. Like ode23s, this solver may be more efficient than ode15s at crude tolerances.

Answer 5

I think this code could work. Try this.

for i =1:n
    t(i +1)=t(i )+dt;
    y(i+1)=solve('y(i+1)=y(i)+dt*f(t(i+1),y(i+1)');
    end 

Answer 6

The code is fine. Just you have to add another loop within the for loop. To check the level of consistency.

if abs((y(i+1) - ynew)/ynew) > 0.0000000001 
    ynew = y(i+1);
    y(i+1) = y(i)+h*f(x(i+1),ynew);
end

I checked for a dummy function and the results were promising.