Help with badly working loop

Hi, hello everyone, I've been working on a method that implements the Backwartd Euler method with the Newton method, amd I have the following code:

function [t,w] = backeuler(f, dfdy, a, b, alpha, N, maxiter, tol)
h = (b-a)/N;
t = a:h:b;
w = t*0;
w(1) = alpha;
converged = false;
for i = 1:N
    w0=w(i);
    wj=w0;
    for j=1:maxiter
        wj=wj-(wj - w0 - h*f(t(i+1),wj)) / (1 - h*dfdy(t(i+1),wj));
        error=(wj - w0 - h*f(t(i+1),wj)) / (1 - h*dfdy(t(i+1),wj));
        fprintf('%d %g\n', j, abs(error));
        if abs(error)<=tol
            converged = true;
            break;
        end
    end
    if converged; break; end
end
fprintf('\n');
if ~converged
    error('No Newton convergence.');
end 

Now, when I terminate:

f = @(t,y) y^2 * (1-y);
dfdy = @(t,y) 2*y - 3*y^2;
a = 0; b = 2000; alpha = 0.9; maxiter = 20; tol = 1e-12;
[t3,w3] = backeuler(f, dfdy, a, b, alpha, 1, maxiter, tol); 
[t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol); 
[t5,w5] = backeuler(f, dfdy, a, b, alpha, 10, maxiter, tol);

I've been getting:

%1st solution
1 0.0270214
2 0.00149354
3 4.46627e-06
4 3.98803e-11
5 7.88716e-18

%2nd solution
1 0.0268375
2 0.00147092
3 4.32538e-06
4 3.7348e-11
5 4.40305e-17

%3rd soltuion
1 0.0266097
2 0.00144316
3 4.15576e-06
4 3.44115e-11
5 1.74336e-17

Which are the three solutions for the givens. But I should be getting more solutions, for example, I only have one solution for [t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol), but I should have four.

%1st solution for [t3,w3] = backeuler(f, dfdy, a, b, alpha, 1, maxiter, tol); 
1 0.128469
2 0.0270214
3 0.00149354
4 4.46627e-006
5 3.98803e-011
6 7.88716e-018

%2nd(1) solution for [t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol); 
1 0.128063
2 0.0268375
3 0.00147092
4 4.32538e-006
5 3.7348e-011
6 4.40305e-017

%2nd(2) soution for [t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol); 
1 0.000249002
2 1.23664e-007
3 3.05022e-014

%2nd(3) solution for [t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol); 
1 6.20646e-007
2 7.68531e-013

%2nd(4) solution for [t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol); 
1 1.54774e-009
2 1.16283e-017

%2nd soltuoin(5)
1 3.8597e-012
2 9.69023e-018

Any help would be appreciated.

[Jan]

You define "error" and try to use the function error() later. This must fail. Avoid shadowing built-in functions by variables. Do you get a corresponding error message?

[Lavorizia Vaughn]

No, I don’t get any error messages.

[Voss]

Seems like you will only ever get one solution because as soon as converged is true, the code breaks out of both for loops.

[Lavorizia Vaughn]

Hi Benjamin, can you please offer a suggestion to my code. What should I change? What should I do? Been stuck on this issue for a couple days now.

[Voss]

I don't know about Backward Euler method, but there are three things I would check on in your code:

1) Check that w is defined correctly. To me, it doesn't seem likely to be correct that the first element of w is alpha and all remaining elements are 0, especially since the elements of w are respectively used as w0.

2) Try swapping the order of the error calculation line and the previous line. And use the error variable in the calculation of wj. That is to say, try changing this:

wj=wj-(wj - w0 - h*f(t(i+1),wj)) / (1 - h*dfdy(t(i+1),wj));

error=(wj - w0 - h*f(t(i+1),wj)) / (1 - h*dfdy(t(i+1),wj));

To this:

error=(wj - w0 - h*f(t(i+1),wj)) / (1 - h*dfdy(t(i+1),wj));

wj=wj-error;

In other words, calculate the error first, then update wj. This change may fix the problem of your error values occuring one iteration sooner than the 'correct' error values you are looking to get.

3) Try moving the converged=false inside the i for-loop, at the top of the loop. This would make the code try to find one solution for each i, so up to N solutions. I don't know if this is right, because I don't know what N is intended for.

(Sorry, code formatting doesn't seem to work when typing on my phone.)

[Lavorizia Vaughn]

Ok I will attempt to make these changes now, and update you shortly.

[Lavorizia Vaughn]

So my code did get better in that the approximations are exactly what they should be, but I am still only getting one solution for the code below, when i should be getting 5.

[t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol);

[Voss]

Regarding my suggestion 3, you'd also need to remove the second break command, i.e., don't break out of the outer loop, but instead continue to find the solution for the next i.

[Lavorizia Vaughn]

N is defined as the the numer of steps or time steps. h depends on N, and h is defined in t=a:h:b;

[Lavorizia Vaughn]

ok I wil make the edit fopr suggestion 3 and update you shortly

[Lavorizia Vaughn]

Yes, im still having the sa,me trouble. Here's my code updated:

function [t,w] = backeuler(f, dfdy, a, b, alpha, N, maxiter, tol)
h = (b-a)/N;
t = a:h:b;
w = t*0;
w(1) = alpha;
for i = 1:N
    converged = false;
    w0=w(i);
    wj=w0;
    for j=1:maxiter
        error=(wj - w0 - h*f(t(i+1),wj)) / (1 - h*dfdy(t(i+1),wj));
        wj=wj-error;
        fprintf('%d %g\n', j, abs(error));
        if abs(error)<=tol
            converged = true;
            break;
        end
    end
    if converged
    end
end
fprintf('\n');
if ~converged
    error('No Newton convergence.');
end 

[Lavorizia Vaughn]

And its returning:

>> [t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol);
1 0.128063
2 0.0268375
3 0.00147092
4 4.32538e-06
5 3.7348e-11
6 4.40305e-17
1 0
1 0
1 0
1 0

When it should be returning:

>> [t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol);
1 0.128063
2 0.0268375
3 0.00147092
4 4.32538e-006
5 3.7348e-011
6 4.40305e-017

1 0.000249002
2 1.23664e-007
3 3.05022e-014

1 6.20646e-007
2 7.68531e-013

1 1.54774e-009
2 1.16283e-017

1 3.8597e-012
2 9.69023e-01

[Voss]

I think you need to use the value wj that has been converged on as the initial wj for the next i-iteration. So, right before or after setting converged = true, add this line: w(i+1) = wj;

[Lavorizia Vaughn]

Dude, that totally worked. You're a godsend, thank you so very much.