Help with badly working loop

7 次查看（过去 30 天）

Lavorizia Vaughn 2021-12-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1601830-help-with-badly-working-loop

评论： Lavorizia Vaughn 2021-12-4

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
0.0270214
0.00149354
4.46627e-06
3.98803e-11
7.88716e-18
%2nd solution
0.0268375
0.00147092
4.32538e-06
3.7348e-11
4.40305e-17
%3rd soltuion
0.0266097
0.00144316
4.15576e-06
3.44115e-11
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); 
0.128469
0.0270214
0.00149354
4.46627e-006
3.98803e-011
7.88716e-018
%2nd(1) solution for [t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol); 
0.128063
0.0268375
0.00147092
4.32538e-006
3.7348e-011
4.40305e-017
%2nd(2) soution for [t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol); 
0.000249002
1.23664e-007
3.05022e-014
%2nd(3) solution for [t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol); 
6.20646e-007
7.68531e-013
%2nd(4) solution for [t4,w4] = backeuler(f, dfdy, a, b, alpha, 5, maxiter, tol); 
1.54774e-009
1.16283e-017
%2nd soltuoin(5)
3.8597e-012
9.69023e-018

Any help would be appreciated.

14 个评论
显示 12更早的评论隐藏 12更早的评论

Voss 2021-12-4

编辑：Voss 2021-12-4

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

Voss 2021-12-4

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 2021-12-4