Mixing
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
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
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
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
asked Dec 24 '18 at 8:50
RebellosRebellos
14.5k31246
asked Dec 24 '18 at 8:50
RebellosRebellos
14.5k31246
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
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 ofstep.
$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 istolnr = 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
});
}
});
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
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
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
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 ofstep.
$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 istolnr = 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
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 ofstep.
$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 istolnr = 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
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
answered Dec 24 '18 at 9:39
LutzLLutzL
57.3k42054
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 ofstep.
$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 istolnr = 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 ofstep.
$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 istolnr = 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
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
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
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
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
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
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
