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

It's hard to try solving an error, if you don't tell us what error or problem you have. — planetmaker, Nov 23 '15 at 13:01
It's recommended that you put your code in the question, so it's still visible if e.g. GitHub goes down. — Serpens, Nov 23 '15 at 13:15

Lutz Lehmann · Accepted Answer · 2015-11-29T15:15:05.993

0

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)
```

edited Nov 29 '15 at 15:15

answered Nov 26 '15 at 12:23

Lutz Lehmann

25,219
2
22
51

Could you help me with that? Exactly with secant.m . I don't understand how to write it for G. – Tiel Wir Nov 29 '15 at 11:44
I rewrote my code with your advice but it still not work. – Tiel Wir Nov 29 '15 at 12:54
There is so much wrong in the general philosophy,... I add a section about that. – Lutz Lehmann Nov 29 '15 at 14:59
Thank you, I rewrote my code with and it is worked now. I still have one question . : Is there is solution for system of equations? – Tiel Wir Dec 05 '15 at 10:28
Yes, there are Newton-like resp. quasi-Newton methods for systems, like the Broyden and adjoint Broyden method that extend the idea of the secant method to systems. Or one uses the exact or approximate Jacobian for Newton's method. – Lutz Lehmann Dec 07 '15 at 08:42

Implicit Euler for stiff equation

1 Answers1