Implicit Euler for stiff equation

Question

Problem: solve stiff different equation
Method: implicit Euler
Plan: I calculate next 'y' by solvin non-linear equation use secant mehod. My function is dy/dx = sin(x+y)

There is right solution . I used newton method

main.m

h=0.01;
x(1)=0;
y_expl(1)=0;  
y_impl(1)=0+h;  
dy(1)=0;
eps=1.0e-6;

for i=1:1000

    x(i+1)=x(i)+h;

    y_impl(i+1)=newton(x(i),y_impl(i),y_impl(i));

    y_expl(i+1)=y_expl(i)+h*f(x(i),y_expl(i));

end

plot(x,y_impl,'r',x,y_expl,'b')
legend('Implicit Euler','Explicit Euler');

newton.m

function [ yn ] = newton( x,y,yi )

    eps=1.0e-6;
    err=1;
    step=0;
    step_max=100;
    h=0.01;
    xn=x+h;

    while (err > eps) && (step < step_max)

        step=step+1;

        yn=y-(F(xn,y,yi,h))/(J(xn,y,h));

        err=abs(y-yn)/(abs(yn)+1.0e-10);

        y=yn;

    end

end

f.m

function [ res ] = f( x,y )

res = sin(x+y);

end

G.m

function [ res ] = J( xn,y,h )

res = h*f(xn,y)-1;

end

F.m

function [ res ] = F( a,y,yn,h )

res = h*f(a,y)-y+yn;

end

Thank for attention

Answer 1

The problem is that you should not be solving F(x,y)=0 but the equation resulting from the implicit Euler step y=y0+h*F(x,y). Thus define

function [res] = G(x,y,y0,h)
    res = y - y0 - h*F(x,y)
end

and use the Newton or secant method for G.

---

General code critique:

• There is no x(0) in matlab.
• implicit Euler is a one-step method, no need to initialize for indices 2 and 3.
• The iteration for the x values is x(i+1)=x(i)+h

---

In the secant method: Known values are x0=x(i), y0=y(i) and h.

• Needed before starting the loop are x1=x0+h and an initial value for y1. This can be taken as the result of an explicit Euler step, y1=y0+h*F(x0,y0) as predictor. The secant method serves as corrector.

• It makes the code more readable if the values of G are computed separately. Note that in G(x1,y,y0,h) the variable is y and the others are fixed parameters. Thus compute G0=G(x1,y0,y0,h) and G1=G(x1,y1,y0,h) for the secant formula y2=y1-G1*(y1-y0)/(G1-G0) or more symmetric y2=(y0*G1-y1*G0)/(G1-G0).

• In principle you could use a generic secant method with interface secant(func, a, b, tol) by calling

  x(i+1) = x(i)+h;
  y(i+1) = secant(@(y) G(x(i+1),y,y(i),h), y(i), y(i)+h*F(x(i),y(i)), delt)