Mixing
How to rewrite Maple code to Matlab, functions fsolve() and solve()?
I have an ordinary differential equation (ODE) -- Van Der Pol oscillator Problem:
y''-2a (1-y^2)y'+y=0,
y(0)=0, y'(0)=0.5, x in [0,10], a =0.025.
The ODE was solved on Maple Software with the fsolve() function.
I need to rewrite the code on Matlab (version R2013a).
My attempt is below.
n = 0
h = 0.1
syms z0; % z in point 0
syms y0; % z in point 0
f0 = 2 * 0.025 * z0 - 2 * 0.025 * ((y0)^2) * z0 - y0
syms z1; % z in point 1
syms y1; % z in point 1
f1 = 2 * 0.025 * z1 - 2 * 0.025 * ((y1)^2) * z1 - y1
syms z2; % z in point 2
syms y2; % z in point 2
f2 = 2 * 0.025 * z2 - 2 * 0.025 * ((y2)^2) * z2 - y2
syms z32; % z in point 3/2
syms y32; % z in point 3/2
f32 = 2 * 0.025 * z32 - 2 * 0.025 * ((y32)^2) * z32 - y32
syms z3; % z in point 3
syms y3; % z in point 3
syms p;
f3 = 2 * p * (1-(y3)^2) * z3 - y3
syms z4; % z in point 4
syms y4; % z in point 4
f4 = 2 * p * (1-(y4)^2) * z4 - y4
syms z72; % z in point 7/2
syms y72; % z in point 7/2
f72 = 2 * p * (1-(y3)^2) * z72 - y72
[c1,c2,c3,c4]=solve(y0, z0, y2 - (-y0 + 2 * y1 + (1/12) * h^2 * f0 + 5/6 * h^2 * f1 + (1/12) * h^2 * f2), y32 - (-1/2 * y0 + 3/2 * y1 + 1/24 *h^2 *f0 + 13/32 * h^2 *f1 - 5/48 * h^2 * f32 + 1/32 * h^2 * f2), -360 * h * z0 - (89 * h^2 * f0 + 186 * h^2 * f1 + 33 * h^2 * f2 - 128 * h^2 * f32 + 360 * y0 -360 * y1))
I have warnings:
Warning: 4 equations in 1 variables.
> In C:Program FilesMATLABR2013atoolboxsymbolicsymbolicsymengine.p>symengine at 56
In mupadengine.mupadengine>mupadengine.evalin at 97
In mupadengine.mupadengine>mupadengine.feval at 150
In solve at 170
Warning: Explicit solution could not be found.
> In solve at 179
Question. Is it possible to solve the problem in Matlab R2013a? Can anyone please give an idea how to rewrite the code?
matlab solver ode maple
asked Dec 24 '18 at 6:59
NickNick
288112
add a comment |
I have an ordinary differential equation (ODE) -- Van Der Pol oscillator Problem:
y''-2a (1-y^2)y'+y=0,
y(0)=0, y'(0)=0.5, x in [0,10], a =0.025.
The ODE was solved on Maple Software with the fsolve() function.
I need to rewrite the code on Matlab (version R2013a).
My attempt is below.
n = 0
h = 0.1
syms z0; % z in point 0
syms y0; % z in point 0
f0 = 2 * 0.025 * z0 - 2 * 0.025 * ((y0)^2) * z0 - y0
syms z1; % z in point 1
syms y1; % z in point 1
f1 = 2 * 0.025 * z1 - 2 * 0.025 * ((y1)^2) * z1 - y1
syms z2; % z in point 2
syms y2; % z in point 2
f2 = 2 * 0.025 * z2 - 2 * 0.025 * ((y2)^2) * z2 - y2
syms z32; % z in point 3/2
syms y32; % z in point 3/2
f32 = 2 * 0.025 * z32 - 2 * 0.025 * ((y32)^2) * z32 - y32
syms z3; % z in point 3
syms y3; % z in point 3
syms p;
f3 = 2 * p * (1-(y3)^2) * z3 - y3
syms z4; % z in point 4
syms y4; % z in point 4
f4 = 2 * p * (1-(y4)^2) * z4 - y4
syms z72; % z in point 7/2
syms y72; % z in point 7/2
f72 = 2 * p * (1-(y3)^2) * z72 - y72
[c1,c2,c3,c4]=solve(y0, z0, y2 - (-y0 + 2 * y1 + (1/12) * h^2 * f0 + 5/6 * h^2 * f1 + (1/12) * h^2 * f2), y32 - (-1/2 * y0 + 3/2 * y1 + 1/24 *h^2 *f0 + 13/32 * h^2 *f1 - 5/48 * h^2 * f32 + 1/32 * h^2 * f2), -360 * h * z0 - (89 * h^2 * f0 + 186 * h^2 * f1 + 33 * h^2 * f2 - 128 * h^2 * f32 + 360 * y0 -360 * y1))
I have warnings:
Warning: 4 equations in 1 variables.
> In C:Program FilesMATLABR2013atoolboxsymbolicsymbolicsymengine.p>symengine at 56
In mupadengine.mupadengine>mupadengine.evalin at 97
In mupadengine.mupadengine>mupadengine.feval at 150
In solve at 170
Warning: Explicit solution could not be found.
> In solve at 179
Question. Is it possible to solve the problem in Matlab R2013a? Can anyone please give an idea how to rewrite the code?
matlab solver ode maple
asked Dec 24 '18 at 6:59
NickNick
288112
What exactly is the method? There has to be some kind of step description that motivates the formulas. Starting withy[n], y[n+1]the aim is to computey[n+2]? This is done by solving an implicit non-linear system for the valuesy[n+3/2], y[n+2], z[n], z[n+1], z[n+3/2], z[n+2]? This is a strange mix of one-step and multi-step approaches. Is there somewhere documentation for this method and its advantages?
– LutzL
Dec 24 '18 at 10:22
add a comment |
I have an ordinary differential equation (ODE) -- Van Der Pol oscillator Problem:
y''-2a (1-y^2)y'+y=0,
y(0)=0, y'(0)=0.5, x in [0,10], a =0.025.
The ODE was solved on Maple Software with the fsolve() function.
I need to rewrite the code on Matlab (version R2013a).
My attempt is below.
n = 0
h = 0.1
syms z0; % z in point 0
syms y0; % z in point 0
f0 = 2 * 0.025 * z0 - 2 * 0.025 * ((y0)^2) * z0 - y0
syms z1; % z in point 1
syms y1; % z in point 1
f1 = 2 * 0.025 * z1 - 2 * 0.025 * ((y1)^2) * z1 - y1
syms z2; % z in point 2
syms y2; % z in point 2
f2 = 2 * 0.025 * z2 - 2 * 0.025 * ((y2)^2) * z2 - y2
syms z32; % z in point 3/2
syms y32; % z in point 3/2
f32 = 2 * 0.025 * z32 - 2 * 0.025 * ((y32)^2) * z32 - y32
syms z3; % z in point 3
syms y3; % z in point 3
syms p;
f3 = 2 * p * (1-(y3)^2) * z3 - y3
syms z4; % z in point 4
syms y4; % z in point 4
f4 = 2 * p * (1-(y4)^2) * z4 - y4
syms z72; % z in point 7/2
syms y72; % z in point 7/2
f72 = 2 * p * (1-(y3)^2) * z72 - y72
[c1,c2,c3,c4]=solve(y0, z0, y2 - (-y0 + 2 * y1 + (1/12) * h^2 * f0 + 5/6 * h^2 * f1 + (1/12) * h^2 * f2), y32 - (-1/2 * y0 + 3/2 * y1 + 1/24 *h^2 *f0 + 13/32 * h^2 *f1 - 5/48 * h^2 * f32 + 1/32 * h^2 * f2), -360 * h * z0 - (89 * h^2 * f0 + 186 * h^2 * f1 + 33 * h^2 * f2 - 128 * h^2 * f32 + 360 * y0 -360 * y1))
I have warnings:
Warning: 4 equations in 1 variables.
> In C:Program FilesMATLABR2013atoolboxsymbolicsymbolicsymengine.p>symengine at 56
In mupadengine.mupadengine>mupadengine.evalin at 97
In mupadengine.mupadengine>mupadengine.feval at 150
In solve at 170
Warning: Explicit solution could not be found.
> In solve at 179
Question. Is it possible to solve the problem in Matlab R2013a? Can anyone please give an idea how to rewrite the code?
matlab solver ode maple
asked Dec 24 '18 at 6:59
NickNick
288112
I have an ordinary differential equation (ODE) -- Van Der Pol oscillator Problem:
y''-2a (1-y^2)y'+y=0,
y(0)=0, y'(0)=0.5, x in [0,10], a =0.025.
The ODE was solved on Maple Software with the fsolve() function.
I need to rewrite the code on Matlab (version R2013a).
My attempt is below.
n = 0
h = 0.1
syms z0; % z in point 0
syms y0; % z in point 0
f0 = 2 * 0.025 * z0 - 2 * 0.025 * ((y0)^2) * z0 - y0
syms z1; % z in point 1
syms y1; % z in point 1
f1 = 2 * 0.025 * z1 - 2 * 0.025 * ((y1)^2) * z1 - y1
syms z2; % z in point 2
syms y2; % z in point 2
f2 = 2 * 0.025 * z2 - 2 * 0.025 * ((y2)^2) * z2 - y2
syms z32; % z in point 3/2
syms y32; % z in point 3/2
f32 = 2 * 0.025 * z32 - 2 * 0.025 * ((y32)^2) * z32 - y32
syms z3; % z in point 3
syms y3; % z in point 3
syms p;
f3 = 2 * p * (1-(y3)^2) * z3 - y3
syms z4; % z in point 4
syms y4; % z in point 4
f4 = 2 * p * (1-(y4)^2) * z4 - y4
syms z72; % z in point 7/2
syms y72; % z in point 7/2
f72 = 2 * p * (1-(y3)^2) * z72 - y72
[c1,c2,c3,c4]=solve(y0, z0, y2 - (-y0 + 2 * y1 + (1/12) * h^2 * f0 + 5/6 * h^2 * f1 + (1/12) * h^2 * f2), y32 - (-1/2 * y0 + 3/2 * y1 + 1/24 *h^2 *f0 + 13/32 * h^2 *f1 - 5/48 * h^2 * f32 + 1/32 * h^2 * f2), -360 * h * z0 - (89 * h^2 * f0 + 186 * h^2 * f1 + 33 * h^2 * f2 - 128 * h^2 * f32 + 360 * y0 -360 * y1))
I have warnings:
Warning: 4 equations in 1 variables.
> In C:Program FilesMATLABR2013atoolboxsymbolicsymbolicsymengine.p>symengine at 56
In mupadengine.mupadengine>mupadengine.evalin at 97
In mupadengine.mupadengine>mupadengine.feval at 150
In solve at 170
Warning: Explicit solution could not be found.
> In solve at 179
Question. Is it possible to solve the problem in Matlab R2013a? Can anyone please give an idea how to rewrite the code?
matlab solver ode maple
matlab solver ode maple
asked Dec 24 '18 at 6:59
NickNick
288112
asked Dec 24 '18 at 6:59
NickNick
288112
asked Dec 24 '18 at 6:59
NickNick
288112
asked Dec 24 '18 at 6:59
NickNick
288112
asked Dec 24 '18 at 6:59
NickNick
288112
288112
What exactly is the method? There has to be some kind of step description that motivates the formulas. Starting withy[n], y[n+1]the aim is to computey[n+2]? This is done by solving an implicit non-linear system for the valuesy[n+3/2], y[n+2], z[n], z[n+1], z[n+3/2], z[n+2]? This is a strange mix of one-step and multi-step approaches. Is there somewhere documentation for this method and its advantages?
– LutzL
Dec 24 '18 at 10:22
add a comment |
What exactly is the method? There has to be some kind of step description that motivates the formulas. Starting withy[n], y[n+1]the aim is to computey[n+2]? This is done by solving an implicit non-linear system for the valuesy[n+3/2], y[n+2], z[n], z[n+1], z[n+3/2], z[n+2]? This is a strange mix of one-step and multi-step approaches. Is there somewhere documentation for this method and its advantages?
– LutzL
Dec 24 '18 at 10:22
What exactly is the method? There has to be some kind of step description that motivates the formulas. Starting with
y[n], y[n+1] the aim is to compute y[n+2]? This is done by solving an implicit non-linear system for the values y[n+3/2], y[n+2], z[n], z[n+1], z[n+3/2], z[n+2]? This is a strange mix of one-step and multi-step approaches. Is there somewhere documentation for this method and its advantages?– LutzL
Dec 24 '18 at 10:22
What exactly is the method? There has to be some kind of step description that motivates the formulas. Starting with
y[n], y[n+1] the aim is to compute y[n+2]? This is done by solving an implicit non-linear system for the values y[n+3/2], y[n+2], z[n], z[n+1], z[n+3/2], z[n+2]? This is a strange mix of one-step and multi-step approaches. Is there somewhere documentation for this method and its advantages?– LutzL
Dec 24 '18 at 10:22
add a comment |
1 Answer
1
active
oldest
votes
The way to find a numerical solution in matlab is to define
VDP_ODE = @(t,u) [ u(2), 2*a*u(2)*(1-u(1)^2) ]; % u = [ y, z ]
(or define it using the full function declaration) and call this with a numerical solver
a=0.025; h=0.1;
t = 0:h:10;
y0 = 0; z0 = 0.5;
u0 = [y0, z0 ];
u = ode45(@VDP_ODE, t, u0)
figure(1); plot(t,u(:,1));
figure(2); plot(u(:,1),u(:,2));
If you want to build your own solver, you should first understand the method. It seems that it is a 4th order method similar but not exactly equal to a linear multi-step method, like the Beeman-Schofield methods introduced to compute molecular dynamics.
The input to each step are the values y(n), y(n+1), the output reduces to y(n+2) with no other values taken over to the next step. Inside the method step additional values are computed for y(n+3/2) and derivatives z at all sample times t(n), t(n+1), t(n+3/2), t(n+2). The aim is to get floating point results, thus it makes no sense to define the equations symbolically, thus adapt the system to the interface of fsolve.
The system that is solved in the Maple code can be simplified by observing that there are 6 equations, thus 6 unknown parameters. The method is based on a 5th degree polynomial with 6 coefficients as unknowns which interpolates the values, first and second derivatives of the solution. This gives 2 equations for the given values y0,y1 and 4 equations for the exact satisfaction of the differential equation at the interpolation points.
Y = [ y0, step0(y0,z0) ]
for n=1:N-2
Y(n+2,:)=stepN(Y(n,:), Y(n+1,:), h)
end
and for the step of the method do something like
%%%% ======== General step ========
%% fit a degree 5 polynomial to values, 1st and 2nd derivatives
%% at t0, t1, t1h=t1+0.5*h, t2
%% p(s) = sum c_k s^k/k!, c0=y1, y(t1+s) ~ p(s)
%% eqn y0 = p(-h)
%% def z0 = p'(-h)
%% eqn f(y0,z0) = p''(-h)
%% def z1 = p'(0)
%% eqn f(y1,z1) = p''(0)
%% eqn y2 = p(h), etc
%%
%% state vector for non-lin solver is u = [y2, c1, c2, ...c5 ]
function y2 = stepN(y0,y1,h)
zz = (y1-y0)/h;
predict = [ y1+h*zz, zz, f(y1,zz), 0, 0, 0];
options = optimset("TolX", 1e-1*h^6, "TolFun", 1e-2*h^6);
u = fsolve(@(u) systemN(u,[y0,y1], h), predict, options);
y2 = u(1);
end
which set some predictor values (only O(h) exact), sets the tolerances for the solver so that the solver error is hopefully smaller than the discretization error which is O(h^6), calls the solver end extracts the wanted value from the solution array.
For the system to solve the function first extracts the constants and variables into named variables for better readability, computes the function values at the given point and then returns the array of the defects in the discretized equations.
function eqn = systemN(u,init,h)
[ y0, y1 ] = num2cell(init){:};
[ y2, c1, c2, c3, c4, c5 ] = num2cell(u){:};
p = @(h) y1 + h*(c1+h/2*(c2+h/3*(c3+h/4*(c4+h/5*c5))));
dp = @(h) c1+h*(c2+h/2*(c3+h/3*(c4+h/4*c5)));
d2p = @(h) c2+h*(c3+h/2*(c4+h/3*c5));
z0 = dp(-h); z1 = c1;
y1h = p(0.5*h); z1h = dp(0.5*h);
z2 = dp(h);
eqn = [ y2-p(h),
y0-p(-h),
f(y0,z0)-d2p(-h),
f(y1,z1)-c2,
f(y1h,z1h)-d2p(0.5*h),
f(y2,z2)-d2p(h) ];
end
answered Dec 24 '18 at 9:18
LutzLLutzL
14.1k21426
add a comment |
Your Answer
StackExchange.ifUsing("editor", function () {
StackExchange.using("externalEditor", function () {
StackExchange.using("snippets", function () {
StackExchange.snippets.init();
});
});
}, "code-snippets");
StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "1"
};
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
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});
}
});
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%2fstackoverflow.com%2fquestions%2f53910201%2fhow-to-rewrite-maple-code-to-matlab-functions-fsolve-and-solve%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
The way to find a numerical solution in matlab is to define
VDP_ODE = @(t,u) [ u(2), 2*a*u(2)*(1-u(1)^2) ]; % u = [ y, z ]
(or define it using the full function declaration) and call this with a numerical solver
a=0.025; h=0.1;
t = 0:h:10;
y0 = 0; z0 = 0.5;
u0 = [y0, z0 ];
u = ode45(@VDP_ODE, t, u0)
figure(1); plot(t,u(:,1));
figure(2); plot(u(:,1),u(:,2));
If you want to build your own solver, you should first understand the method. It seems that it is a 4th order method similar but not exactly equal to a linear multi-step method, like the Beeman-Schofield methods introduced to compute molecular dynamics.
The input to each step are the values y(n), y(n+1), the output reduces to y(n+2) with no other values taken over to the next step. Inside the method step additional values are computed for y(n+3/2) and derivatives z at all sample times t(n), t(n+1), t(n+3/2), t(n+2). The aim is to get floating point results, thus it makes no sense to define the equations symbolically, thus adapt the system to the interface of fsolve.
The system that is solved in the Maple code can be simplified by observing that there are 6 equations, thus 6 unknown parameters. The method is based on a 5th degree polynomial with 6 coefficients as unknowns which interpolates the values, first and second derivatives of the solution. This gives 2 equations for the given values y0,y1 and 4 equations for the exact satisfaction of the differential equation at the interpolation points.
Y = [ y0, step0(y0,z0) ]
for n=1:N-2
Y(n+2,:)=stepN(Y(n,:), Y(n+1,:), h)
end
and for the step of the method do something like
%%%% ======== General step ========
%% fit a degree 5 polynomial to values, 1st and 2nd derivatives
%% at t0, t1, t1h=t1+0.5*h, t2
%% p(s) = sum c_k s^k/k!, c0=y1, y(t1+s) ~ p(s)
%% eqn y0 = p(-h)
%% def z0 = p'(-h)
%% eqn f(y0,z0) = p''(-h)
%% def z1 = p'(0)
%% eqn f(y1,z1) = p''(0)
%% eqn y2 = p(h), etc
%%
%% state vector for non-lin solver is u = [y2, c1, c2, ...c5 ]
function y2 = stepN(y0,y1,h)
zz = (y1-y0)/h;
predict = [ y1+h*zz, zz, f(y1,zz), 0, 0, 0];
options = optimset("TolX", 1e-1*h^6, "TolFun", 1e-2*h^6);
u = fsolve(@(u) systemN(u,[y0,y1], h), predict, options);
y2 = u(1);
end
which set some predictor values (only O(h) exact), sets the tolerances for the solver so that the solver error is hopefully smaller than the discretization error which is O(h^6), calls the solver end extracts the wanted value from the solution array.
For the system to solve the function first extracts the constants and variables into named variables for better readability, computes the function values at the given point and then returns the array of the defects in the discretized equations.
function eqn = systemN(u,init,h)
[ y0, y1 ] = num2cell(init){:};
[ y2, c1, c2, c3, c4, c5 ] = num2cell(u){:};
p = @(h) y1 + h*(c1+h/2*(c2+h/3*(c3+h/4*(c4+h/5*c5))));
dp = @(h) c1+h*(c2+h/2*(c3+h/3*(c4+h/4*c5)));
d2p = @(h) c2+h*(c3+h/2*(c4+h/3*c5));
z0 = dp(-h); z1 = c1;
y1h = p(0.5*h); z1h = dp(0.5*h);
z2 = dp(h);
eqn = [ y2-p(h),
y0-p(-h),
f(y0,z0)-d2p(-h),
f(y1,z1)-c2,
f(y1h,z1h)-d2p(0.5*h),
f(y2,z2)-d2p(h) ];
end
answered Dec 24 '18 at 9:18
LutzLLutzL
14.1k21426
add a comment |
The way to find a numerical solution in matlab is to define
VDP_ODE = @(t,u) [ u(2), 2*a*u(2)*(1-u(1)^2) ]; % u = [ y, z ]
(or define it using the full function declaration) and call this with a numerical solver
a=0.025; h=0.1;
t = 0:h:10;
y0 = 0; z0 = 0.5;
u0 = [y0, z0 ];
u = ode45(@VDP_ODE, t, u0)
figure(1); plot(t,u(:,1));
figure(2); plot(u(:,1),u(:,2));
If you want to build your own solver, you should first understand the method. It seems that it is a 4th order method similar but not exactly equal to a linear multi-step method, like the Beeman-Schofield methods introduced to compute molecular dynamics.
The input to each step are the values y(n), y(n+1), the output reduces to y(n+2) with no other values taken over to the next step. Inside the method step additional values are computed for y(n+3/2) and derivatives z at all sample times t(n), t(n+1), t(n+3/2), t(n+2). The aim is to get floating point results, thus it makes no sense to define the equations symbolically, thus adapt the system to the interface of fsolve.
The system that is solved in the Maple code can be simplified by observing that there are 6 equations, thus 6 unknown parameters. The method is based on a 5th degree polynomial with 6 coefficients as unknowns which interpolates the values, first and second derivatives of the solution. This gives 2 equations for the given values y0,y1 and 4 equations for the exact satisfaction of the differential equation at the interpolation points.
Y = [ y0, step0(y0,z0) ]
for n=1:N-2
Y(n+2,:)=stepN(Y(n,:), Y(n+1,:), h)
end
and for the step of the method do something like
%%%% ======== General step ========
%% fit a degree 5 polynomial to values, 1st and 2nd derivatives
%% at t0, t1, t1h=t1+0.5*h, t2
%% p(s) = sum c_k s^k/k!, c0=y1, y(t1+s) ~ p(s)
%% eqn y0 = p(-h)
%% def z0 = p'(-h)
%% eqn f(y0,z0) = p''(-h)
%% def z1 = p'(0)
%% eqn f(y1,z1) = p''(0)
%% eqn y2 = p(h), etc
%%
%% state vector for non-lin solver is u = [y2, c1, c2, ...c5 ]
function y2 = stepN(y0,y1,h)
zz = (y1-y0)/h;
predict = [ y1+h*zz, zz, f(y1,zz), 0, 0, 0];
options = optimset("TolX", 1e-1*h^6, "TolFun", 1e-2*h^6);
u = fsolve(@(u) systemN(u,[y0,y1], h), predict, options);
y2 = u(1);
end
which set some predictor values (only O(h) exact), sets the tolerances for the solver so that the solver error is hopefully smaller than the discretization error which is O(h^6), calls the solver end extracts the wanted value from the solution array.
For the system to solve the function first extracts the constants and variables into named variables for better readability, computes the function values at the given point and then returns the array of the defects in the discretized equations.
function eqn = systemN(u,init,h)
[ y0, y1 ] = num2cell(init){:};
[ y2, c1, c2, c3, c4, c5 ] = num2cell(u){:};
p = @(h) y1 + h*(c1+h/2*(c2+h/3*(c3+h/4*(c4+h/5*c5))));
dp = @(h) c1+h*(c2+h/2*(c3+h/3*(c4+h/4*c5)));
d2p = @(h) c2+h*(c3+h/2*(c4+h/3*c5));
z0 = dp(-h); z1 = c1;
y1h = p(0.5*h); z1h = dp(0.5*h);
z2 = dp(h);
eqn = [ y2-p(h),
y0-p(-h),
f(y0,z0)-d2p(-h),
f(y1,z1)-c2,
f(y1h,z1h)-d2p(0.5*h),
f(y2,z2)-d2p(h) ];
end
answered Dec 24 '18 at 9:18
LutzLLutzL
14.1k21426
add a comment |
The way to find a numerical solution in matlab is to define
VDP_ODE = @(t,u) [ u(2), 2*a*u(2)*(1-u(1)^2) ]; % u = [ y, z ]
(or define it using the full function declaration) and call this with a numerical solver
a=0.025; h=0.1;
t = 0:h:10;
y0 = 0; z0 = 0.5;
u0 = [y0, z0 ];
u = ode45(@VDP_ODE, t, u0)
figure(1); plot(t,u(:,1));
figure(2); plot(u(:,1),u(:,2));
If you want to build your own solver, you should first understand the method. It seems that it is a 4th order method similar but not exactly equal to a linear multi-step method, like the Beeman-Schofield methods introduced to compute molecular dynamics.
The input to each step are the values y(n), y(n+1), the output reduces to y(n+2) with no other values taken over to the next step. Inside the method step additional values are computed for y(n+3/2) and derivatives z at all sample times t(n), t(n+1), t(n+3/2), t(n+2). The aim is to get floating point results, thus it makes no sense to define the equations symbolically, thus adapt the system to the interface of fsolve.
The system that is solved in the Maple code can be simplified by observing that there are 6 equations, thus 6 unknown parameters. The method is based on a 5th degree polynomial with 6 coefficients as unknowns which interpolates the values, first and second derivatives of the solution. This gives 2 equations for the given values y0,y1 and 4 equations for the exact satisfaction of the differential equation at the interpolation points.
Y = [ y0, step0(y0,z0) ]
for n=1:N-2
Y(n+2,:)=stepN(Y(n,:), Y(n+1,:), h)
end
and for the step of the method do something like
%%%% ======== General step ========
%% fit a degree 5 polynomial to values, 1st and 2nd derivatives
%% at t0, t1, t1h=t1+0.5*h, t2
%% p(s) = sum c_k s^k/k!, c0=y1, y(t1+s) ~ p(s)
%% eqn y0 = p(-h)
%% def z0 = p'(-h)
%% eqn f(y0,z0) = p''(-h)
%% def z1 = p'(0)
%% eqn f(y1,z1) = p''(0)
%% eqn y2 = p(h), etc
%%
%% state vector for non-lin solver is u = [y2, c1, c2, ...c5 ]
function y2 = stepN(y0,y1,h)
zz = (y1-y0)/h;
predict = [ y1+h*zz, zz, f(y1,zz), 0, 0, 0];
options = optimset("TolX", 1e-1*h^6, "TolFun", 1e-2*h^6);
u = fsolve(@(u) systemN(u,[y0,y1], h), predict, options);
y2 = u(1);
end
which set some predictor values (only O(h) exact), sets the tolerances for the solver so that the solver error is hopefully smaller than the discretization error which is O(h^6), calls the solver end extracts the wanted value from the solution array.
For the system to solve the function first extracts the constants and variables into named variables for better readability, computes the function values at the given point and then returns the array of the defects in the discretized equations.
function eqn = systemN(u,init,h)
[ y0, y1 ] = num2cell(init){:};
[ y2, c1, c2, c3, c4, c5 ] = num2cell(u){:};
p = @(h) y1 + h*(c1+h/2*(c2+h/3*(c3+h/4*(c4+h/5*c5))));
dp = @(h) c1+h*(c2+h/2*(c3+h/3*(c4+h/4*c5)));
d2p = @(h) c2+h*(c3+h/2*(c4+h/3*c5));
z0 = dp(-h); z1 = c1;
y1h = p(0.5*h); z1h = dp(0.5*h);
z2 = dp(h);
eqn = [ y2-p(h),
y0-p(-h),
f(y0,z0)-d2p(-h),
f(y1,z1)-c2,
f(y1h,z1h)-d2p(0.5*h),
f(y2,z2)-d2p(h) ];
end
answered Dec 24 '18 at 9:18
LutzLLutzL
14.1k21426
The way to find a numerical solution in matlab is to define
VDP_ODE = @(t,u) [ u(2), 2*a*u(2)*(1-u(1)^2) ]; % u = [ y, z ]
(or define it using the full function declaration) and call this with a numerical solver
a=0.025; h=0.1;
t = 0:h:10;
y0 = 0; z0 = 0.5;
u0 = [y0, z0 ];
u = ode45(@VDP_ODE, t, u0)
figure(1); plot(t,u(:,1));
figure(2); plot(u(:,1),u(:,2));
If you want to build your own solver, you should first understand the method. It seems that it is a 4th order method similar but not exactly equal to a linear multi-step method, like the Beeman-Schofield methods introduced to compute molecular dynamics.
The input to each step are the values y(n), y(n+1), the output reduces to y(n+2) with no other values taken over to the next step. Inside the method step additional values are computed for y(n+3/2) and derivatives z at all sample times t(n), t(n+1), t(n+3/2), t(n+2). The aim is to get floating point results, thus it makes no sense to define the equations symbolically, thus adapt the system to the interface of fsolve.
The system that is solved in the Maple code can be simplified by observing that there are 6 equations, thus 6 unknown parameters. The method is based on a 5th degree polynomial with 6 coefficients as unknowns which interpolates the values, first and second derivatives of the solution. This gives 2 equations for the given values y0,y1 and 4 equations for the exact satisfaction of the differential equation at the interpolation points.
Y = [ y0, step0(y0,z0) ]
for n=1:N-2
Y(n+2,:)=stepN(Y(n,:), Y(n+1,:), h)
end
and for the step of the method do something like
%%%% ======== General step ========
%% fit a degree 5 polynomial to values, 1st and 2nd derivatives
%% at t0, t1, t1h=t1+0.5*h, t2
%% p(s) = sum c_k s^k/k!, c0=y1, y(t1+s) ~ p(s)
%% eqn y0 = p(-h)
%% def z0 = p'(-h)
%% eqn f(y0,z0) = p''(-h)
%% def z1 = p'(0)
%% eqn f(y1,z1) = p''(0)
%% eqn y2 = p(h), etc
%%
%% state vector for non-lin solver is u = [y2, c1, c2, ...c5 ]
function y2 = stepN(y0,y1,h)
zz = (y1-y0)/h;
predict = [ y1+h*zz, zz, f(y1,zz), 0, 0, 0];
options = optimset("TolX", 1e-1*h^6, "TolFun", 1e-2*h^6);
u = fsolve(@(u) systemN(u,[y0,y1], h), predict, options);
y2 = u(1);
end
which set some predictor values (only O(h) exact), sets the tolerances for the solver so that the solver error is hopefully smaller than the discretization error which is O(h^6), calls the solver end extracts the wanted value from the solution array.
For the system to solve the function first extracts the constants and variables into named variables for better readability, computes the function values at the given point and then returns the array of the defects in the discretized equations.
function eqn = systemN(u,init,h)
[ y0, y1 ] = num2cell(init){:};
[ y2, c1, c2, c3, c4, c5 ] = num2cell(u){:};
p = @(h) y1 + h*(c1+h/2*(c2+h/3*(c3+h/4*(c4+h/5*c5))));
dp = @(h) c1+h*(c2+h/2*(c3+h/3*(c4+h/4*c5)));
d2p = @(h) c2+h*(c3+h/2*(c4+h/3*c5));
z0 = dp(-h); z1 = c1;
y1h = p(0.5*h); z1h = dp(0.5*h);
z2 = dp(h);
eqn = [ y2-p(h),
y0-p(-h),
f(y0,z0)-d2p(-h),
f(y1,z1)-c2,
f(y1h,z1h)-d2p(0.5*h),
f(y2,z2)-d2p(h) ];
end
answered Dec 24 '18 at 9:18
LutzLLutzL
14.1k21426
edited Jan 1 at 18:13
answered Dec 24 '18 at 9:18
LutzLLutzL
14.1k21426
answered Dec 24 '18 at 9:18
LutzLLutzL
14.1k21426
answered Dec 24 '18 at 9:18
LutzLLutzL
14.1k21426
14.1k21426
add a comment |
add a comment |
Thanks for contributing an answer to Stack Overflow!
- 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.
To learn more, see our tips on writing great answers.
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%2fstackoverflow.com%2fquestions%2f53910201%2fhow-to-rewrite-maple-code-to-matlab-functions-fsolve-and-solve%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
What exactly is the method? There has to be some kind of step description that motivates the formulas. Starting with
y[n], y[n+1]the aim is to computey[n+2]? This is done by solving an implicit non-linear system for the valuesy[n+3/2], y[n+2], z[n], z[n+1], z[n+3/2], z[n+2]? This is a strange mix of one-step and multi-step approaches. Is there somewhere documentation for this method and its advantages?– LutzL
Dec 24 '18 at 10:22