Matlab code for $y_{n+2} - y_n = hleft[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_nright]$

1

$begingroup$

I am looking for a Matlab code for the $2$ stage multistep method :
$$y_{n+2} - y_n = hleft[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_nright]$$
Note that $f_n = f(x_n,y_n)$. Now, since $f_{n+2} = f(x_{n+2},y_{n+2})$ and $y_{n+2}$ is also the output of the method, cleary this is an implicit method. That means that Newton-Raphson shall be used for the calculation of the final output. The point of that code is to use it in a paper-work I am assembling for some conclusions over methods. Theoritically, I have proven that this is a method of order $4$, meaning that : $|T_n| leq ch^4$. I want to prove this experimentally by using a certain method, but this can't be done if I do not have the code for the given method.

Since I am not that good in programming, I am kindly requesting some help regarding the code. I would really appreciate any elaborated algorithm. Thanks in advance !

My attempt at a code using the rk4 function already defined for the first steps :

function [tout, yout] = newcorrect(FunFcn,t0,tfinal,step,y0)

maxiter = 1000;

tolnr   = 1e-9;

diffdelta = 1e-6;

stages=2;

[tout,yout]=rk4(FunFcn,t0,t0+(stages-1)*step,step,y0);

tout=tout(1:stages);

yout=yout(1:stages);

t = tout(stages);

y = yout(stages).';

% The main loop

 while abs(t- tfinal)> 1e-6 

    if t + step > tfinal, step = tfinal - t; end

    t = t + step;

    yn0 = y;

    ynf = yn0;

    yn = inf;

    iter = 0;

    while (abs(yn - ynf)>= tolnr) && (iter < maxiter)

        df = 1/diffdelta * (feval(FunFcn,t, yn0+diffdelta) - feval(FunFcn, t, yn0));

        yn = yn0 - 1/(1/3*step*df - 1) * (4/3*step*feval(FunFcn,tout(end),yout(end)) + 1/3*step*feval(FunFcn,tout(end-1),yout(end-1)) + 1/3*step*feval(FunFcn, t, yn0) -yn0 + yout(end-1));

        ynf = yn0;

        yn0 = yn;

        iter = iter + 1;

    end

    y = yn;

    tout = [tout; t];

    yout = [yout; y.'];

 end

end

ordinary-differential-equations numerical-methods matlab newton-raphson

share|cite|improve this question

edited Dec 24 '18 at 18:10

Rebellos

asked Dec 24 '18 at 8:50

RebellosRebellos

14.5k31246

$endgroup$

add a comment | 

1

$begingroup$

I am looking for a Matlab code for the $2$ stage multistep method :
$$y_{n+2} - y_n = hleft[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_nright]$$
Note that $f_n = f(x_n,y_n)$. Now, since $f_{n+2} = f(x_{n+2},y_{n+2})$ and $y_{n+2}$ is also the output of the method, cleary this is an implicit method. That means that Newton-Raphson shall be used for the calculation of the final output. The point of that code is to use it in a paper-work I am assembling for some conclusions over methods. Theoritically, I have proven that this is a method of order $4$, meaning that : $|T_n| leq ch^4$. I want to prove this experimentally by using a certain method, but this can't be done if I do not have the code for the given method.

Since I am not that good in programming, I am kindly requesting some help regarding the code. I would really appreciate any elaborated algorithm. Thanks in advance !

My attempt at a code using the rk4 function already defined for the first steps :

function [tout, yout] = newcorrect(FunFcn,t0,tfinal,step,y0)

maxiter = 1000;

tolnr   = 1e-9;

diffdelta = 1e-6;

stages=2;

[tout,yout]=rk4(FunFcn,t0,t0+(stages-1)*step,step,y0);

tout=tout(1:stages);

yout=yout(1:stages);

t = tout(stages);

y = yout(stages).';

% The main loop

 while abs(t- tfinal)> 1e-6 

    if t + step > tfinal, step = tfinal - t; end

    t = t + step;

    yn0 = y;

    ynf = yn0;

    yn = inf;

    iter = 0;

    while (abs(yn - ynf)>= tolnr) && (iter < maxiter)

        df = 1/diffdelta * (feval(FunFcn,t, yn0+diffdelta) - feval(FunFcn, t, yn0));

        yn = yn0 - 1/(1/3*step*df - 1) * (4/3*step*feval(FunFcn,tout(end),yout(end)) + 1/3*step*feval(FunFcn,tout(end-1),yout(end-1)) + 1/3*step*feval(FunFcn, t, yn0) -yn0 + yout(end-1));

        ynf = yn0;

        yn0 = yn;

        iter = iter + 1;

    end

    y = yn;

    tout = [tout; t];

    yout = [yout; y.'];

 end

end

ordinary-differential-equations numerical-methods matlab newton-raphson

share|cite|improve this question

edited Dec 24 '18 at 18:10

Rebellos

asked Dec 24 '18 at 8:50

RebellosRebellos

14.5k31246

$endgroup$

add a comment | 

1

1

1

$begingroup$

I am looking for a Matlab code for the $2$ stage multistep method :
$$y_{n+2} - y_n = hleft[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_nright]$$
Note that $f_n = f(x_n,y_n)$. Now, since $f_{n+2} = f(x_{n+2},y_{n+2})$ and $y_{n+2}$ is also the output of the method, cleary this is an implicit method. That means that Newton-Raphson shall be used for the calculation of the final output. The point of that code is to use it in a paper-work I am assembling for some conclusions over methods. Theoritically, I have proven that this is a method of order $4$, meaning that : $|T_n| leq ch^4$. I want to prove this experimentally by using a certain method, but this can't be done if I do not have the code for the given method.

Since I am not that good in programming, I am kindly requesting some help regarding the code. I would really appreciate any elaborated algorithm. Thanks in advance !

My attempt at a code using the rk4 function already defined for the first steps :

function [tout, yout] = newcorrect(FunFcn,t0,tfinal,step,y0)

maxiter = 1000;

tolnr   = 1e-9;

diffdelta = 1e-6;

stages=2;

[tout,yout]=rk4(FunFcn,t0,t0+(stages-1)*step,step,y0);

tout=tout(1:stages);

yout=yout(1:stages);

t = tout(stages);

y = yout(stages).';

% The main loop

 while abs(t- tfinal)> 1e-6 

    if t + step > tfinal, step = tfinal - t; end

    t = t + step;

    yn0 = y;

    ynf = yn0;

    yn = inf;

    iter = 0;

    while (abs(yn - ynf)>= tolnr) && (iter < maxiter)

        df = 1/diffdelta * (feval(FunFcn,t, yn0+diffdelta) - feval(FunFcn, t, yn0));

        yn = yn0 - 1/(1/3*step*df - 1) * (4/3*step*feval(FunFcn,tout(end),yout(end)) + 1/3*step*feval(FunFcn,tout(end-1),yout(end-1)) + 1/3*step*feval(FunFcn, t, yn0) -yn0 + yout(end-1));

        ynf = yn0;

        yn0 = yn;

        iter = iter + 1;

    end

    y = yn;

    tout = [tout; t];

    yout = [yout; y.'];

 end

end

ordinary-differential-equations numerical-methods matlab newton-raphson

share|cite|improve this question

edited Dec 24 '18 at 18:10

Rebellos

asked Dec 24 '18 at 8:50

RebellosRebellos

14.5k31246

$endgroup$

I am looking for a Matlab code for the $2$ stage multistep method :
$$y_{n+2} - y_n = hleft[(1/3)f_{n+2} + (4/3)f_{n+1} + (1/3)f_nright]$$
Note that $f_n = f(x_n,y_n)$. Now, since $f_{n+2} = f(x_{n+2},y_{n+2})$ and $y_{n+2}$ is also the output of the method, cleary this is an implicit method. That means that Newton-Raphson shall be used for the calculation of the final output. The point of that code is to use it in a paper-work I am assembling for some conclusions over methods. Theoritically, I have proven that this is a method of order $4$, meaning that : $|T_n| leq ch^4$. I want to prove this experimentally by using a certain method, but this can't be done if I do not have the code for the given method.

Since I am not that good in programming, I am kindly requesting some help regarding the code. I would really appreciate any elaborated algorithm. Thanks in advance !

My attempt at a code using the rk4 function already defined for the first steps :

function [tout, yout] = newcorrect(FunFcn,t0,tfinal,step,y0)

maxiter = 1000;

tolnr   = 1e-9;

diffdelta = 1e-6;

stages=2;

[tout,yout]=rk4(FunFcn,t0,t0+(stages-1)*step,step,y0);

tout=tout(1:stages);

yout=yout(1:stages);

t = tout(stages);

y = yout(stages).';

% The main loop

 while abs(t- tfinal)> 1e-6 

    if t + step > tfinal, step = tfinal - t; end

    t = t + step;

    yn0 = y;

    ynf = yn0;

    yn = inf;

    iter = 0;

    while (abs(yn - ynf)>= tolnr) && (iter < maxiter)

        df = 1/diffdelta * (feval(FunFcn,t, yn0+diffdelta) - feval(FunFcn, t, yn0));

        yn = yn0 - 1/(1/3*step*df - 1) * (4/3*step*feval(FunFcn,tout(end),yout(end)) + 1/3*step*feval(FunFcn,tout(end-1),yout(end-1)) + 1/3*step*feval(FunFcn, t, yn0) -yn0 + yout(end-1));

        ynf = yn0;

        yn0 = yn;

        iter = iter + 1;

    end

    y = yn;

    tout = [tout; t];

    yout = [yout; y.'];

 end

end

ordinary-differential-equations numerical-methods matlab newton-raphson

ordinary-differential-equations numerical-methods matlab newton-raphson

share|cite|improve this question

edited Dec 24 '18 at 18:10

Rebellos

asked Dec 24 '18 at 8:50

RebellosRebellos

14.5k31246

share|cite|improve this question

edited Dec 24 '18 at 18:10

Rebellos

asked Dec 24 '18 at 8:50

RebellosRebellos

14.5k31246

share|cite|improve this question

share|cite|improve this question

edited Dec 24 '18 at 18:10

Rebellos

edited Dec 24 '18 at 18:10

Rebellos

edited Dec 24 '18 at 18:10

Rebellos

asked Dec 24 '18 at 8:50

RebellosRebellos

14.5k31246

asked Dec 24 '18 at 8:50

RebellosRebellos

14.5k31246

asked Dec 24 '18 at 8:50

RebellosRebellos

14.5k31246

14.5k31246

add a comment | 

add a comment | 

1 Answer
1

active

oldest

votes

3

$begingroup$

Using fsolve to short-cut solving the implicit equation without explicitly coding a Newton method gives the following implementation

function [T, Y] = multistep(FunFcn,t0,tfinal,step,y0)

  % initialize arrays

  N = 1;

  T = [ t0 ];

  Y = [ y0 ];

  F = [ FunFcn(t0, y0) ];



  % perform an RK3 step to have initial error O(h^4)

  % to get the requisite two samples to start the multistep method,

  % or use RK4 to stay with the O(h^5) local error

  N = N+1;

  T(N) = T(N-1)+step;

  Y(N,:) = RK3step(FunFcn,T(N-1),step,Y(N-1,:));

  F(N,:) = FunFcn(T(N), Y(N,:));



  % set error tolerances for implicit solver, 

  % the implicit error (times step size) should 

  % contribute less than the discretization error

  options = optimset("TolX", 1e-3*step^4, "TolFun", 1e-6*step^4);



  % main loop

  while abs(T(N)+0.6*step - tfinal) > 0.501*step 

    N = N+1;

    T(N) = T(N-1)+step ;

    % perform the multi-step method

    % y2 -y0 = h/3*(f0+4*f1+f2)

    % y2 - h/3*f(t2,y2) = b = y0 + h/3*(f0+4*f1)

    b = Y(N-2,:) + step/3*(F(N-2,:)+4*F(N-1,:));

    % solve the implicit equation using the built-in fsolve

    % Coding your own solver would allow a better management 

    % of the Jacobian approximation

    correct = @(y) ( y-step/3*FunFcn(T(N), y) - b );

    % predict = Y(N-2,:)+2*step*F(N-1,:);

    % use a 4th order predictor

    predict = 5*Y(N-2,:)-4*Y(N-1,:)+2*step*(2*F(N-1,:)+F(N-2,:));

    Y(N,:) = fsolve(correct, predict, options );

    F(N,:) = FunFcn(T(N), Y(N,:));

  end



  % Now tfinal is between T(N) + 0.1*step and T(N)+1.1*step

  % compute the final value with an RK3 step to add no more than an O(h^4) error

  N = N+1;

  T(N) = tfinal;

  Y(N,:) = RK3step(FunFcn,T(N-1),T(N)-T(N-1), Y(N-1,:));

end    

Using RK3 to handle the first and last steps is the minimum required, it is probably better to use RK4 to have the error contributions negligible against the global error order 4.

The Runge-Kutta step functions implement the well-known methods:

function y1=RK3step(FunFcn, t0, step, y0)

  k1 = step*FunFcn(t0,y0);

  k2 = step*FunFcn(t0+0.5*step, y0+0.5*k1);

  k3 = step*FunFcn(t0+    step, y0+2*k2-k1);

  y1 = y0 + (k1+4*k2+k3)/6;

end



function y1=RK4step(FunFcn, t0, step, y0)

  k1 = step*FunFcn(t0         , y0       );

  k2 = step*FunFcn(t0+0.5*step, y0+0.5*k1);

  k3 = step*FunFcn(t0+0.5*step, y0+0.5*k2);

  k4 = step*FunFcn(t0+    step, y0+    k3);

  y1 = y0 + (k1+2*k2+2*k3+k4)/6;

end

Now test the error order on a problem with known exact solution

%% test problem y''+sin(y) = r = p''+sin(p), y(0)=p(0), y'(0)=p'(0)

%% p=sin(t), p(0)=0, p'(0)=1, r(t)=-sin(t)+sin(sin(t))



TestODE = @(t,u) [ u(2), -sin(u(1))-sin(t)+sin(sin(t)) ];



t0 = 0; tf = 2;

[ T, Y ] = multistep(TestODE, t0,tf,0.2, [ 0, 1 ]);



printf("values found y1(1)=%.15f, y2(1)=%15fn", Y(end,1), Y(end,2));

printf("exact values y1(1)=%.15f, y2(1)=%15fn", sin(tf)  , cos(tf)  );

figure(1);

plot(Y(:,1), Y(:,2));

And to confirm the order, logarithmically plot the error against the step size

H = 0.1.^[1:0.5:3];

E1=; E2=;

for h = H

  [ T, Y ] = multistep(@(t,u)TestODE(t,u),t0,tf,h,[ 0, 1 ]);

  E1(end+1) = max(abs( Y(:,1)'-sin(T) ));

  E2(end+1) = max(abs( Y(:,2)'-cos(T) ));

end

figure(2)

loglog(H,E1,H,E2)

end

share|cite|improve this answer

edited Jan 5 at 13:04

answered Dec 24 '18 at 9:39

LutzLLutzL

57.3k42054

$endgroup$

• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Hi and thanks a lot for your elaborative input ! I faced some difficulties inputing functions. I added my own attempt at a code, which despite working and providing very good approximations for solutions, it seems to be of order $1$ and not $4$ as it should be, which troubles me a lot. Would you mind taking a look
  $endgroup$
  – Rebellos
  Dec 24 '18 at 17:49
  
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  It might be that the last step is the problem here. You could compute the error leaving out the last values. The method is only exact when the interval length is an integer multiple of step.
  $endgroup$
  – LutzL
  Dec 24 '18 at 18:31
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Appreciate the input to start off. I am solving the ODE : $y' = frac{t^2+y^2}{2ty}$ over the interval $[1,3]$ with $y(1) = 2$. Then, by the directions of the exercise, I should make a step $h = 2^{-k}$ and the fraction of the maximum consecutive absolute errors for each $h$should tend to $2^{-p}$ where $p$ is the order of the method, thus $1/16$. Using my code, I all of the values of the fractions yielded indeed tend to $1/16$ BUT the last one for $h=2^{-9}$ and $h=2^{-10}$ all of a sudden is $0.24$. How is that explained ? Or is there something wrong in my code?
  $endgroup$
  – Rebellos
  Dec 24 '18 at 18:35
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Your error tolerance for the Newton method is tolnr = 1e-9, the global error in general can not be smaller than that. Reduce that proportional to $h^4$ to get an adapted error tolerance. // It should be sufficient to compute the derivative approximation once per integration step. The simplified Newton method has then a contraction quotient of $O(h^2)$. If you use the Euler step as predictor, the Newton iterates will then have errors $O(h^2), ~ O(h^4), ~ O(h^6)$ and the last one is in general accurate enough.
  $endgroup$
  – LutzL
  Dec 24 '18 at 19:08
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Great elaboration, thanks ! I can safely say I wasn't aware of these small details. I tried lowering the tolerances in general and that last term became lower, so I see what you meant in practice as well ! A TON of thanks for all the dedication to the post and the great answer ( I will work over your code and implement it as well ). Merry Christmas !
  $endgroup$
  – Rebellos
  Dec 24 '18 at 19:47
  

  
  
  
  

 | 
show 1 more comment

Your Answer

StackExchange.ifUsing("editor", function () {
return StackExchange.using("mathjaxEditing", function () {
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
});
});
}, "mathjax-editing");

StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "69"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});

function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
noCode: true, onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});


}
});

Thanks for contributing an answer to Mathematics Stack Exchange!

• Please be sure to answer the question. Provide details and share your research!

But avoid …

• Asking for help, clarification, or responding to other answers.


• Making statements based on opinion; back them up with references or personal experience.

Use MathJax to format equations. MathJax reference.

To learn more, see our tips on writing great answers.

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name

Email

Required, but never shown

StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f3051078%2fmatlab-code-for-y-n2-y-n-h-left1-3f-n2-4-3f-n1-1-3f-n-r%23new-answer', 'question_page');
}
);

Post as a guest

Name

Email

Required, but never shown

1 Answer
1

active

oldest

votes

1 Answer
1

active

oldest

votes

active

oldest

votes

active

oldest

votes

3

$begingroup$

Using fsolve to short-cut solving the implicit equation without explicitly coding a Newton method gives the following implementation

function [T, Y] = multistep(FunFcn,t0,tfinal,step,y0)

  % initialize arrays

  N = 1;

  T = [ t0 ];

  Y = [ y0 ];

  F = [ FunFcn(t0, y0) ];



  % perform an RK3 step to have initial error O(h^4)

  % to get the requisite two samples to start the multistep method,

  % or use RK4 to stay with the O(h^5) local error

  N = N+1;

  T(N) = T(N-1)+step;

  Y(N,:) = RK3step(FunFcn,T(N-1),step,Y(N-1,:));

  F(N,:) = FunFcn(T(N), Y(N,:));



  % set error tolerances for implicit solver, 

  % the implicit error (times step size) should 

  % contribute less than the discretization error

  options = optimset("TolX", 1e-3*step^4, "TolFun", 1e-6*step^4);



  % main loop

  while abs(T(N)+0.6*step - tfinal) > 0.501*step 

    N = N+1;

    T(N) = T(N-1)+step ;

    % perform the multi-step method

    % y2 -y0 = h/3*(f0+4*f1+f2)

    % y2 - h/3*f(t2,y2) = b = y0 + h/3*(f0+4*f1)

    b = Y(N-2,:) + step/3*(F(N-2,:)+4*F(N-1,:));

    % solve the implicit equation using the built-in fsolve

    % Coding your own solver would allow a better management 

    % of the Jacobian approximation

    correct = @(y) ( y-step/3*FunFcn(T(N), y) - b );

    % predict = Y(N-2,:)+2*step*F(N-1,:);

    % use a 4th order predictor

    predict = 5*Y(N-2,:)-4*Y(N-1,:)+2*step*(2*F(N-1,:)+F(N-2,:));

    Y(N,:) = fsolve(correct, predict, options );

    F(N,:) = FunFcn(T(N), Y(N,:));

  end



  % Now tfinal is between T(N) + 0.1*step and T(N)+1.1*step

  % compute the final value with an RK3 step to add no more than an O(h^4) error

  N = N+1;

  T(N) = tfinal;

  Y(N,:) = RK3step(FunFcn,T(N-1),T(N)-T(N-1), Y(N-1,:));

end    

Using RK3 to handle the first and last steps is the minimum required, it is probably better to use RK4 to have the error contributions negligible against the global error order 4.

The Runge-Kutta step functions implement the well-known methods:

function y1=RK3step(FunFcn, t0, step, y0)

  k1 = step*FunFcn(t0,y0);

  k2 = step*FunFcn(t0+0.5*step, y0+0.5*k1);

  k3 = step*FunFcn(t0+    step, y0+2*k2-k1);

  y1 = y0 + (k1+4*k2+k3)/6;

end



function y1=RK4step(FunFcn, t0, step, y0)

  k1 = step*FunFcn(t0         , y0       );

  k2 = step*FunFcn(t0+0.5*step, y0+0.5*k1);

  k3 = step*FunFcn(t0+0.5*step, y0+0.5*k2);

  k4 = step*FunFcn(t0+    step, y0+    k3);

  y1 = y0 + (k1+2*k2+2*k3+k4)/6;

end

Now test the error order on a problem with known exact solution

%% test problem y''+sin(y) = r = p''+sin(p), y(0)=p(0), y'(0)=p'(0)

%% p=sin(t), p(0)=0, p'(0)=1, r(t)=-sin(t)+sin(sin(t))



TestODE = @(t,u) [ u(2), -sin(u(1))-sin(t)+sin(sin(t)) ];



t0 = 0; tf = 2;

[ T, Y ] = multistep(TestODE, t0,tf,0.2, [ 0, 1 ]);



printf("values found y1(1)=%.15f, y2(1)=%15fn", Y(end,1), Y(end,2));

printf("exact values y1(1)=%.15f, y2(1)=%15fn", sin(tf)  , cos(tf)  );

figure(1);

plot(Y(:,1), Y(:,2));

And to confirm the order, logarithmically plot the error against the step size

H = 0.1.^[1:0.5:3];

E1=; E2=;

for h = H

  [ T, Y ] = multistep(@(t,u)TestODE(t,u),t0,tf,h,[ 0, 1 ]);

  E1(end+1) = max(abs( Y(:,1)'-sin(T) ));

  E2(end+1) = max(abs( Y(:,2)'-cos(T) ));

end

figure(2)

loglog(H,E1,H,E2)

end

share|cite|improve this answer

edited Jan 5 at 13:04

answered Dec 24 '18 at 9:39

LutzLLutzL

57.3k42054

$endgroup$

• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Hi and thanks a lot for your elaborative input ! I faced some difficulties inputing functions. I added my own attempt at a code, which despite working and providing very good approximations for solutions, it seems to be of order $1$ and not $4$ as it should be, which troubles me a lot. Would you mind taking a look
  $endgroup$
  – Rebellos
  Dec 24 '18 at 17:49
  
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  It might be that the last step is the problem here. You could compute the error leaving out the last values. The method is only exact when the interval length is an integer multiple of step.
  $endgroup$
  – LutzL
  Dec 24 '18 at 18:31
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Appreciate the input to start off. I am solving the ODE : $y' = frac{t^2+y^2}{2ty}$ over the interval $[1,3]$ with $y(1) = 2$. Then, by the directions of the exercise, I should make a step $h = 2^{-k}$ and the fraction of the maximum consecutive absolute errors for each $h$should tend to $2^{-p}$ where $p$ is the order of the method, thus $1/16$. Using my code, I all of the values of the fractions yielded indeed tend to $1/16$ BUT the last one for $h=2^{-9}$ and $h=2^{-10}$ all of a sudden is $0.24$. How is that explained ? Or is there something wrong in my code?
  $endgroup$
  – Rebellos
  Dec 24 '18 at 18:35
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Your error tolerance for the Newton method is tolnr = 1e-9, the global error in general can not be smaller than that. Reduce that proportional to $h^4$ to get an adapted error tolerance. // It should be sufficient to compute the derivative approximation once per integration step. The simplified Newton method has then a contraction quotient of $O(h^2)$. If you use the Euler step as predictor, the Newton iterates will then have errors $O(h^2), ~ O(h^4), ~ O(h^6)$ and the last one is in general accurate enough.
  $endgroup$
  – LutzL
  Dec 24 '18 at 19:08
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Great elaboration, thanks ! I can safely say I wasn't aware of these small details. I tried lowering the tolerances in general and that last term became lower, so I see what you meant in practice as well ! A TON of thanks for all the dedication to the post and the great answer ( I will work over your code and implement it as well ). Merry Christmas !
  $endgroup$
  – Rebellos
  Dec 24 '18 at 19:47
  

  
  
  
  

 | 
show 1 more comment

3

$begingroup$

Using fsolve to short-cut solving the implicit equation without explicitly coding a Newton method gives the following implementation

function [T, Y] = multistep(FunFcn,t0,tfinal,step,y0)

  % initialize arrays

  N = 1;

  T = [ t0 ];

  Y = [ y0 ];

  F = [ FunFcn(t0, y0) ];



  % perform an RK3 step to have initial error O(h^4)

  % to get the requisite two samples to start the multistep method,

  % or use RK4 to stay with the O(h^5) local error

  N = N+1;

  T(N) = T(N-1)+step;

  Y(N,:) = RK3step(FunFcn,T(N-1),step,Y(N-1,:));

  F(N,:) = FunFcn(T(N), Y(N,:));



  % set error tolerances for implicit solver, 

  % the implicit error (times step size) should 

  % contribute less than the discretization error

  options = optimset("TolX", 1e-3*step^4, "TolFun", 1e-6*step^4);



  % main loop

  while abs(T(N)+0.6*step - tfinal) > 0.501*step 

    N = N+1;

    T(N) = T(N-1)+step ;

    % perform the multi-step method

    % y2 -y0 = h/3*(f0+4*f1+f2)

    % y2 - h/3*f(t2,y2) = b = y0 + h/3*(f0+4*f1)

    b = Y(N-2,:) + step/3*(F(N-2,:)+4*F(N-1,:));

    % solve the implicit equation using the built-in fsolve

    % Coding your own solver would allow a better management 

    % of the Jacobian approximation

    correct = @(y) ( y-step/3*FunFcn(T(N), y) - b );

    % predict = Y(N-2,:)+2*step*F(N-1,:);

    % use a 4th order predictor

    predict = 5*Y(N-2,:)-4*Y(N-1,:)+2*step*(2*F(N-1,:)+F(N-2,:));

    Y(N,:) = fsolve(correct, predict, options );

    F(N,:) = FunFcn(T(N), Y(N,:));

  end



  % Now tfinal is between T(N) + 0.1*step and T(N)+1.1*step

  % compute the final value with an RK3 step to add no more than an O(h^4) error

  N = N+1;

  T(N) = tfinal;

  Y(N,:) = RK3step(FunFcn,T(N-1),T(N)-T(N-1), Y(N-1,:));

end    

Using RK3 to handle the first and last steps is the minimum required, it is probably better to use RK4 to have the error contributions negligible against the global error order 4.

The Runge-Kutta step functions implement the well-known methods:

function y1=RK3step(FunFcn, t0, step, y0)

  k1 = step*FunFcn(t0,y0);

  k2 = step*FunFcn(t0+0.5*step, y0+0.5*k1);

  k3 = step*FunFcn(t0+    step, y0+2*k2-k1);

  y1 = y0 + (k1+4*k2+k3)/6;

end



function y1=RK4step(FunFcn, t0, step, y0)

  k1 = step*FunFcn(t0         , y0       );

  k2 = step*FunFcn(t0+0.5*step, y0+0.5*k1);

  k3 = step*FunFcn(t0+0.5*step, y0+0.5*k2);

  k4 = step*FunFcn(t0+    step, y0+    k3);

  y1 = y0 + (k1+2*k2+2*k3+k4)/6;

end

Now test the error order on a problem with known exact solution

%% test problem y''+sin(y) = r = p''+sin(p), y(0)=p(0), y'(0)=p'(0)

%% p=sin(t), p(0)=0, p'(0)=1, r(t)=-sin(t)+sin(sin(t))



TestODE = @(t,u) [ u(2), -sin(u(1))-sin(t)+sin(sin(t)) ];



t0 = 0; tf = 2;

[ T, Y ] = multistep(TestODE, t0,tf,0.2, [ 0, 1 ]);



printf("values found y1(1)=%.15f, y2(1)=%15fn", Y(end,1), Y(end,2));

printf("exact values y1(1)=%.15f, y2(1)=%15fn", sin(tf)  , cos(tf)  );

figure(1);

plot(Y(:,1), Y(:,2));

And to confirm the order, logarithmically plot the error against the step size

H = 0.1.^[1:0.5:3];

E1=; E2=;

for h = H

  [ T, Y ] = multistep(@(t,u)TestODE(t,u),t0,tf,h,[ 0, 1 ]);

  E1(end+1) = max(abs( Y(:,1)'-sin(T) ));

  E2(end+1) = max(abs( Y(:,2)'-cos(T) ));

end

figure(2)

loglog(H,E1,H,E2)

end

share|cite|improve this answer

edited Jan 5 at 13:04

answered Dec 24 '18 at 9:39

LutzLLutzL

57.3k42054

$endgroup$

• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Hi and thanks a lot for your elaborative input ! I faced some difficulties inputing functions. I added my own attempt at a code, which despite working and providing very good approximations for solutions, it seems to be of order $1$ and not $4$ as it should be, which troubles me a lot. Would you mind taking a look
  $endgroup$
  – Rebellos
  Dec 24 '18 at 17:49
  
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  It might be that the last step is the problem here. You could compute the error leaving out the last values. The method is only exact when the interval length is an integer multiple of step.
  $endgroup$
  – LutzL
  Dec 24 '18 at 18:31
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Appreciate the input to start off. I am solving the ODE : $y' = frac{t^2+y^2}{2ty}$ over the interval $[1,3]$ with $y(1) = 2$. Then, by the directions of the exercise, I should make a step $h = 2^{-k}$ and the fraction of the maximum consecutive absolute errors for each $h$should tend to $2^{-p}$ where $p$ is the order of the method, thus $1/16$. Using my code, I all of the values of the fractions yielded indeed tend to $1/16$ BUT the last one for $h=2^{-9}$ and $h=2^{-10}$ all of a sudden is $0.24$. How is that explained ? Or is there something wrong in my code?
  $endgroup$
  – Rebellos
  Dec 24 '18 at 18:35
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Your error tolerance for the Newton method is tolnr = 1e-9, the global error in general can not be smaller than that. Reduce that proportional to $h^4$ to get an adapted error tolerance. // It should be sufficient to compute the derivative approximation once per integration step. The simplified Newton method has then a contraction quotient of $O(h^2)$. If you use the Euler step as predictor, the Newton iterates will then have errors $O(h^2), ~ O(h^4), ~ O(h^6)$ and the last one is in general accurate enough.
  $endgroup$
  – LutzL
  Dec 24 '18 at 19:08
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Great elaboration, thanks ! I can safely say I wasn't aware of these small details. I tried lowering the tolerances in general and that last term became lower, so I see what you meant in practice as well ! A TON of thanks for all the dedication to the post and the great answer ( I will work over your code and implement it as well ). Merry Christmas !
  $endgroup$
  – Rebellos
  Dec 24 '18 at 19:47
  

  
  
  
  

 | 
show 1 more comment

3

3

3

$begingroup$

Using fsolve to short-cut solving the implicit equation without explicitly coding a Newton method gives the following implementation

function [T, Y] = multistep(FunFcn,t0,tfinal,step,y0)

  % initialize arrays

  N = 1;

  T = [ t0 ];

  Y = [ y0 ];

  F = [ FunFcn(t0, y0) ];



  % perform an RK3 step to have initial error O(h^4)

  % to get the requisite two samples to start the multistep method,

  % or use RK4 to stay with the O(h^5) local error

  N = N+1;

  T(N) = T(N-1)+step;

  Y(N,:) = RK3step(FunFcn,T(N-1),step,Y(N-1,:));

  F(N,:) = FunFcn(T(N), Y(N,:));



  % set error tolerances for implicit solver, 

  % the implicit error (times step size) should 

  % contribute less than the discretization error

  options = optimset("TolX", 1e-3*step^4, "TolFun", 1e-6*step^4);



  % main loop

  while abs(T(N)+0.6*step - tfinal) > 0.501*step 

    N = N+1;

    T(N) = T(N-1)+step ;

    % perform the multi-step method

    % y2 -y0 = h/3*(f0+4*f1+f2)

    % y2 - h/3*f(t2,y2) = b = y0 + h/3*(f0+4*f1)

    b = Y(N-2,:) + step/3*(F(N-2,:)+4*F(N-1,:));

    % solve the implicit equation using the built-in fsolve

    % Coding your own solver would allow a better management 

    % of the Jacobian approximation

    correct = @(y) ( y-step/3*FunFcn(T(N), y) - b );

    % predict = Y(N-2,:)+2*step*F(N-1,:);

    % use a 4th order predictor

    predict = 5*Y(N-2,:)-4*Y(N-1,:)+2*step*(2*F(N-1,:)+F(N-2,:));

    Y(N,:) = fsolve(correct, predict, options );

    F(N,:) = FunFcn(T(N), Y(N,:));

  end



  % Now tfinal is between T(N) + 0.1*step and T(N)+1.1*step

  % compute the final value with an RK3 step to add no more than an O(h^4) error

  N = N+1;

  T(N) = tfinal;

  Y(N,:) = RK3step(FunFcn,T(N-1),T(N)-T(N-1), Y(N-1,:));

end    

Using RK3 to handle the first and last steps is the minimum required, it is probably better to use RK4 to have the error contributions negligible against the global error order 4.

The Runge-Kutta step functions implement the well-known methods:

function y1=RK3step(FunFcn, t0, step, y0)

  k1 = step*FunFcn(t0,y0);

  k2 = step*FunFcn(t0+0.5*step, y0+0.5*k1);

  k3 = step*FunFcn(t0+    step, y0+2*k2-k1);

  y1 = y0 + (k1+4*k2+k3)/6;

end



function y1=RK4step(FunFcn, t0, step, y0)

  k1 = step*FunFcn(t0         , y0       );

  k2 = step*FunFcn(t0+0.5*step, y0+0.5*k1);

  k3 = step*FunFcn(t0+0.5*step, y0+0.5*k2);

  k4 = step*FunFcn(t0+    step, y0+    k3);

  y1 = y0 + (k1+2*k2+2*k3+k4)/6;

end

Now test the error order on a problem with known exact solution

%% test problem y''+sin(y) = r = p''+sin(p), y(0)=p(0), y'(0)=p'(0)

%% p=sin(t), p(0)=0, p'(0)=1, r(t)=-sin(t)+sin(sin(t))



TestODE = @(t,u) [ u(2), -sin(u(1))-sin(t)+sin(sin(t)) ];



t0 = 0; tf = 2;

[ T, Y ] = multistep(TestODE, t0,tf,0.2, [ 0, 1 ]);



printf("values found y1(1)=%.15f, y2(1)=%15fn", Y(end,1), Y(end,2));

printf("exact values y1(1)=%.15f, y2(1)=%15fn", sin(tf)  , cos(tf)  );

figure(1);

plot(Y(:,1), Y(:,2));

And to confirm the order, logarithmically plot the error against the step size

H = 0.1.^[1:0.5:3];

E1=; E2=;

for h = H

  [ T, Y ] = multistep(@(t,u)TestODE(t,u),t0,tf,h,[ 0, 1 ]);

  E1(end+1) = max(abs( Y(:,1)'-sin(T) ));

  E2(end+1) = max(abs( Y(:,2)'-cos(T) ));

end

figure(2)

loglog(H,E1,H,E2)

end

share|cite|improve this answer

edited Jan 5 at 13:04

answered Dec 24 '18 at 9:39

LutzLLutzL

57.3k42054

$endgroup$

Using fsolve to short-cut solving the implicit equation without explicitly coding a Newton method gives the following implementation

function [T, Y] = multistep(FunFcn,t0,tfinal,step,y0)

  % initialize arrays

  N = 1;

  T = [ t0 ];

  Y = [ y0 ];

  F = [ FunFcn(t0, y0) ];



  % perform an RK3 step to have initial error O(h^4)

  % to get the requisite two samples to start the multistep method,

  % or use RK4 to stay with the O(h^5) local error

  N = N+1;

  T(N) = T(N-1)+step;

  Y(N,:) = RK3step(FunFcn,T(N-1),step,Y(N-1,:));

  F(N,:) = FunFcn(T(N), Y(N,:));



  % set error tolerances for implicit solver, 

  % the implicit error (times step size) should 

  % contribute less than the discretization error

  options = optimset("TolX", 1e-3*step^4, "TolFun", 1e-6*step^4);



  % main loop

  while abs(T(N)+0.6*step - tfinal) > 0.501*step 

    N = N+1;

    T(N) = T(N-1)+step ;

    % perform the multi-step method

    % y2 -y0 = h/3*(f0+4*f1+f2)

    % y2 - h/3*f(t2,y2) = b = y0 + h/3*(f0+4*f1)

    b = Y(N-2,:) + step/3*(F(N-2,:)+4*F(N-1,:));

    % solve the implicit equation using the built-in fsolve

    % Coding your own solver would allow a better management 

    % of the Jacobian approximation

    correct = @(y) ( y-step/3*FunFcn(T(N), y) - b );

    % predict = Y(N-2,:)+2*step*F(N-1,:);

    % use a 4th order predictor

    predict = 5*Y(N-2,:)-4*Y(N-1,:)+2*step*(2*F(N-1,:)+F(N-2,:));

    Y(N,:) = fsolve(correct, predict, options );

    F(N,:) = FunFcn(T(N), Y(N,:));

  end



  % Now tfinal is between T(N) + 0.1*step and T(N)+1.1*step

  % compute the final value with an RK3 step to add no more than an O(h^4) error

  N = N+1;

  T(N) = tfinal;

  Y(N,:) = RK3step(FunFcn,T(N-1),T(N)-T(N-1), Y(N-1,:));

end    

Using RK3 to handle the first and last steps is the minimum required, it is probably better to use RK4 to have the error contributions negligible against the global error order 4.

The Runge-Kutta step functions implement the well-known methods:

function y1=RK3step(FunFcn, t0, step, y0)

  k1 = step*FunFcn(t0,y0);

  k2 = step*FunFcn(t0+0.5*step, y0+0.5*k1);

  k3 = step*FunFcn(t0+    step, y0+2*k2-k1);

  y1 = y0 + (k1+4*k2+k3)/6;

end



function y1=RK4step(FunFcn, t0, step, y0)

  k1 = step*FunFcn(t0         , y0       );

  k2 = step*FunFcn(t0+0.5*step, y0+0.5*k1);

  k3 = step*FunFcn(t0+0.5*step, y0+0.5*k2);

  k4 = step*FunFcn(t0+    step, y0+    k3);

  y1 = y0 + (k1+2*k2+2*k3+k4)/6;

end

Now test the error order on a problem with known exact solution

%% test problem y''+sin(y) = r = p''+sin(p), y(0)=p(0), y'(0)=p'(0)

%% p=sin(t), p(0)=0, p'(0)=1, r(t)=-sin(t)+sin(sin(t))



TestODE = @(t,u) [ u(2), -sin(u(1))-sin(t)+sin(sin(t)) ];



t0 = 0; tf = 2;

[ T, Y ] = multistep(TestODE, t0,tf,0.2, [ 0, 1 ]);



printf("values found y1(1)=%.15f, y2(1)=%15fn", Y(end,1), Y(end,2));

printf("exact values y1(1)=%.15f, y2(1)=%15fn", sin(tf)  , cos(tf)  );

figure(1);

plot(Y(:,1), Y(:,2));

And to confirm the order, logarithmically plot the error against the step size

H = 0.1.^[1:0.5:3];

E1=; E2=;

for h = H

  [ T, Y ] = multistep(@(t,u)TestODE(t,u),t0,tf,h,[ 0, 1 ]);

  E1(end+1) = max(abs( Y(:,1)'-sin(T) ));

  E2(end+1) = max(abs( Y(:,2)'-cos(T) ));

end

figure(2)

loglog(H,E1,H,E2)

end

share|cite|improve this answer

edited Jan 5 at 13:04

answered Dec 24 '18 at 9:39

LutzLLutzL

57.3k42054

share|cite|improve this answer

share|cite|improve this answer

edited Jan 5 at 13:04

edited Jan 5 at 13:04

edited Jan 5 at 13:04

answered Dec 24 '18 at 9:39

LutzLLutzL

57.3k42054

answered Dec 24 '18 at 9:39

LutzLLutzL

57.3k42054

answered Dec 24 '18 at 9:39

LutzLLutzL

57.3k42054

57.3k42054

• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Hi and thanks a lot for your elaborative input ! I faced some difficulties inputing functions. I added my own attempt at a code, which despite working and providing very good approximations for solutions, it seems to be of order $1$ and not $4$ as it should be, which troubles me a lot. Would you mind taking a look
  $endgroup$
  – Rebellos
  Dec 24 '18 at 17:49
  
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  It might be that the last step is the problem here. You could compute the error leaving out the last values. The method is only exact when the interval length is an integer multiple of step.
  $endgroup$
  – LutzL
  Dec 24 '18 at 18:31
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Appreciate the input to start off. I am solving the ODE : $y' = frac{t^2+y^2}{2ty}$ over the interval $[1,3]$ with $y(1) = 2$. Then, by the directions of the exercise, I should make a step $h = 2^{-k}$ and the fraction of the maximum consecutive absolute errors for each $h$should tend to $2^{-p}$ where $p$ is the order of the method, thus $1/16$. Using my code, I all of the values of the fractions yielded indeed tend to $1/16$ BUT the last one for $h=2^{-9}$ and $h=2^{-10}$ all of a sudden is $0.24$. How is that explained ? Or is there something wrong in my code?
  $endgroup$
  – Rebellos
  Dec 24 '18 at 18:35
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Your error tolerance for the Newton method is tolnr = 1e-9, the global error in general can not be smaller than that. Reduce that proportional to $h^4$ to get an adapted error tolerance. // It should be sufficient to compute the derivative approximation once per integration step. The simplified Newton method has then a contraction quotient of $O(h^2)$. If you use the Euler step as predictor, the Newton iterates will then have errors $O(h^2), ~ O(h^4), ~ O(h^6)$ and the last one is in general accurate enough.
  $endgroup$
  – LutzL
  Dec 24 '18 at 19:08
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Great elaboration, thanks ! I can safely say I wasn't aware of these small details. I tried lowering the tolerances in general and that last term became lower, so I see what you meant in practice as well ! A TON of thanks for all the dedication to the post and the great answer ( I will work over your code and implement it as well ). Merry Christmas !
  $endgroup$
  – Rebellos
  Dec 24 '18 at 19:47
  

  
  
  
  

 | 
show 1 more comment

• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Hi and thanks a lot for your elaborative input ! I faced some difficulties inputing functions. I added my own attempt at a code, which despite working and providing very good approximations for solutions, it seems to be of order $1$ and not $4$ as it should be, which troubles me a lot. Would you mind taking a look
  $endgroup$
  – Rebellos
  Dec 24 '18 at 17:49
  
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  It might be that the last step is the problem here. You could compute the error leaving out the last values. The method is only exact when the interval length is an integer multiple of step.
  $endgroup$
  – LutzL
  Dec 24 '18 at 18:31
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Appreciate the input to start off. I am solving the ODE : $y' = frac{t^2+y^2}{2ty}$ over the interval $[1,3]$ with $y(1) = 2$. Then, by the directions of the exercise, I should make a step $h = 2^{-k}$ and the fraction of the maximum consecutive absolute errors for each $h$should tend to $2^{-p}$ where $p$ is the order of the method, thus $1/16$. Using my code, I all of the values of the fractions yielded indeed tend to $1/16$ BUT the last one for $h=2^{-9}$ and $h=2^{-10}$ all of a sudden is $0.24$. How is that explained ? Or is there something wrong in my code?
  $endgroup$
  – Rebellos
  Dec 24 '18 at 18:35
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Your error tolerance for the Newton method is tolnr = 1e-9, the global error in general can not be smaller than that. Reduce that proportional to $h^4$ to get an adapted error tolerance. // It should be sufficient to compute the derivative approximation once per integration step. The simplified Newton method has then a contraction quotient of $O(h^2)$. If you use the Euler step as predictor, the Newton iterates will then have errors $O(h^2), ~ O(h^4), ~ O(h^6)$ and the last one is in general accurate enough.
  $endgroup$
  – LutzL
  Dec 24 '18 at 19:08
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  $begingroup$
  Great elaboration, thanks ! I can safely say I wasn't aware of these small details. I tried lowering the tolerances in general and that last term became lower, so I see what you meant in practice as well ! A TON of thanks for all the dedication to the post and the great answer ( I will work over your code and implement it as well ). Merry Christmas !
  $endgroup$
  – Rebellos
  Dec 24 '18 at 19:47
  

  
  
  
  

$begingroup$
Hi and thanks a lot for your elaborative input ! I faced some difficulties inputing functions. I added my own attempt at a code, which despite working and providing very good approximations for solutions, it seems to be of order $1$ and not $4$ as it should be, which troubles me a lot. Would you mind taking a look
$endgroup$
– Rebellos
Dec 24 '18 at 17:49

$begingroup$
Hi and thanks a lot for your elaborative input ! I faced some difficulties inputing functions. I added my own attempt at a code, which despite working and providing very good approximations for solutions, it seems to be of order $1$ and not $4$ as it should be, which troubles me a lot. Would you mind taking a look
$endgroup$
– Rebellos
Dec 24 '18 at 17:49

$begingroup$
It might be that the last step is the problem here. You could compute the error leaving out the last values. The method is only exact when the interval length is an integer multiple of step.
$endgroup$
– LutzL
Dec 24 '18 at 18:31

$begingroup$
It might be that the last step is the problem here. You could compute the error leaving out the last values. The method is only exact when the interval length is an integer multiple of step.
$endgroup$
– LutzL
Dec 24 '18 at 18:31

$begingroup$
Appreciate the input to start off. I am solving the ODE : $y' = frac{t^2+y^2}{2ty}$ over the interval $[1,3]$ with $y(1) = 2$. Then, by the directions of the exercise, I should make a step $h = 2^{-k}$ and the fraction of the maximum consecutive absolute errors for each $h$should tend to $2^{-p}$ where $p$ is the order of the method, thus $1/16$. Using my code, I all of the values of the fractions yielded indeed tend to $1/16$ BUT the last one for $h=2^{-9}$ and $h=2^{-10}$ all of a sudden is $0.24$. How is that explained ? Or is there something wrong in my code?
$endgroup$
– Rebellos
Dec 24 '18 at 18:35

$begingroup$
Appreciate the input to start off. I am solving the ODE : $y' = frac{t^2+y^2}{2ty}$ over the interval $[1,3]$ with $y(1) = 2$. Then, by the directions of the exercise, I should make a step $h = 2^{-k}$ and the fraction of the maximum consecutive absolute errors for each $h$should tend to $2^{-p}$ where $p$ is the order of the method, thus $1/16$. Using my code, I all of the values of the fractions yielded indeed tend to $1/16$ BUT the last one for $h=2^{-9}$ and $h=2^{-10}$ all of a sudden is $0.24$. How is that explained ? Or is there something wrong in my code?
$endgroup$
– Rebellos
Dec 24 '18 at 18:35

$begingroup$
Your error tolerance for the Newton method is tolnr = 1e-9, the global error in general can not be smaller than that. Reduce that proportional to $h^4$ to get an adapted error tolerance. // It should be sufficient to compute the derivative approximation once per integration step. The simplified Newton method has then a contraction quotient of $O(h^2)$. If you use the Euler step as predictor, the Newton iterates will then have errors $O(h^2), ~ O(h^4), ~ O(h^6)$ and the last one is in general accurate enough.
$endgroup$
– LutzL
Dec 24 '18 at 19:08

$begingroup$
Your error tolerance for the Newton method is tolnr = 1e-9, the global error in general can not be smaller than that. Reduce that proportional to $h^4$ to get an adapted error tolerance. // It should be sufficient to compute the derivative approximation once per integration step. The simplified Newton method has then a contraction quotient of $O(h^2)$. If you use the Euler step as predictor, the Newton iterates will then have errors $O(h^2), ~ O(h^4), ~ O(h^6)$ and the last one is in general accurate enough.
$endgroup$
– LutzL
Dec 24 '18 at 19:08

$begingroup$
Great elaboration, thanks ! I can safely say I wasn't aware of these small details. I tried lowering the tolerances in general and that last term became lower, so I see what you meant in practice as well ! A TON of thanks for all the dedication to the post and the great answer ( I will work over your code and implement it as well ). Merry Christmas !
$endgroup$
– Rebellos
Dec 24 '18 at 19:47

$begingroup$
Great elaboration, thanks ! I can safely say I wasn't aware of these small details. I tried lowering the tolerances in general and that last term became lower, so I see what you meant in practice as well ! A TON of thanks for all the dedication to the post and the great answer ( I will work over your code and implement it as well ). Merry Christmas !
$endgroup$
– Rebellos
Dec 24 '18 at 19:47

 | 
show 1 more comment

Thanks for contributing an answer to Mathematics Stack Exchange!

• Please be sure to answer the question. Provide details and share your research!

But avoid …

• Asking for help, clarification, or responding to other answers.


• Making statements based on opinion; back them up with references or personal experience.

Use MathJax to format equations. MathJax reference.

To learn more, see our tips on writing great answers.

draft saved

draft discarded

Thanks for contributing an answer to Mathematics Stack Exchange!

• Please be sure to answer the question. Provide details and share your research!

But avoid …

• Asking for help, clarification, or responding to other answers.


• Making statements based on opinion; back them up with references or personal experience.

Use MathJax to format equations. MathJax reference.

To learn more, see our tips on writing great answers.

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name

Email

Required, but never shown

StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f3051078%2fmatlab-code-for-y-n2-y-n-h-left1-3f-n2-4-3f-n1-1-3f-n-r%23new-answer', 'question_page');
}
);

Post as a guest

Name

Email

Required, but never shown

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name

Email

Required, but never shown

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name

Email

Required, but never shown

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name

Email

Required, but never shown

Name

Name

Email

Required, but never shown

Email

Required, but never shown

Email

Required, but never shown

Email

Required, but never shown

Name

Name

Email

Required, but never shown

Email

Required, but never shown

Email

Required, but never shown

Email

Required, but never shown

This page is only for reference, If you need detailed information, please check here