Ways to speed up user implemented RK4Speed up Numerical IntegrationSpeed of convergence for NIntegrateTough Calculation, novice mathematica userNumerical integration's speedNumerical integral speedImprove the speed of Gaussian quadrature integrationSolving an unstable BVP numerically, accurately and efficientlyHow to speed up integral of results of PDE modelSolve BVP involving user defined functionUser defined ArcTan function
Teaching indefinite integrals that require special-casing
What to do with wrong results in talks?
How do I define a right arrow with bar in LaTeX?
Is a roofing delivery truck likely to crack my driveway slab?
How will losing mobility of one hand affect my career as a programmer?
Why are on-board computers allowed to change controls without notifying the pilots?
What is the intuitive meaning of having a linear relationship between the logs of two variables?
Have I saved too much for retirement so far?
Can somebody explain Brexit in a few child-proof sentences?
Can criminal fraud exist without damages?
voltage of sounds of mp3files
How can I replace every global instance of "x[2]" with "x_2"
Is it correct to write "is not focus on"?
Why does John Bercow say “unlock” after reading out the results of a vote?
Is there any reason not to eat food that's been dropped on the surface of the moon?
Failed to fetch jessie backports repository
How can a jailer prevent the Forge Cleric's Artisan's Blessing from being used?
Student evaluations of teaching assistants
I'm in charge of equipment buying but no one's ever happy with what I choose. How to fix this?
Hide Select Output from T-SQL
Star/Wye electrical connection math symbol
Can a monster with multiattack use this ability if they are missing a limb?
Why did Kant, Hegel, and Adorno leave some words and phrases in the Greek alphabet?
Go Pregnant or Go Home
Ways to speed up user implemented RK4
Speed up Numerical IntegrationSpeed of convergence for NIntegrateTough Calculation, novice mathematica userNumerical integration's speedNumerical integral speedImprove the speed of Gaussian quadrature integrationSolving an unstable BVP numerically, accurately and efficientlyHow to speed up integral of results of PDE modelSolve BVP involving user defined functionUser defined ArcTan function
3
$begingroup$
So, I've implemented RK4, and I'm wondering what I can do to make it more efficient? What I've got so far is below. I wish to still record all steps. I think AppendTo is doing the most damage to the time, is there a faster alternative?
rk4[f_, variables_, valtinit_, tinit_, tfinal_, nsteps_] :=
Module[table, xlist, ylist, step, k1, k2, k3, k4,
xlist = tinit;
step = N[(tfinal - tinit)/(nsteps)];
ylist = valtinit;
table = xlist, ylist;
Table[
k1 = step* f /. MapThread[Rule, variables, ylist]; (*
Equivalent to step* f/.Thread[Rule[variables,ylist]]*)
k2 = step*f /. MapThread[Rule, variables, k1/2 + ylist];
k3 = step*f /. MapThread[Rule, variables, k2/2 + ylist];
k4 = step*f /. MapThread[Rule, variables, k3 + ylist];
ylist += 1/6 (k1 + 2 (k2 + k3) + k4);
xlist += step;
AppendTo[table, xlist, ylist];
xlist, ylist, nsteps];
table
];
Example Input:
funclist = -x + y, x - y;
initials = 1, 2;
variables = x, y;
init = 0;
final = 200;
nstep = 20000;
approx = rk4[funclist, variables, initials, init, final, nstep]//AbsoluteTiming;
3.59932,...
I'd love some suggestions!
differential-equations numerical-integration
asked 4 hours ago
ShinaolordShinaolord
808
$endgroup$
|
show 1 more comment
3
$begingroup$
So, I've implemented RK4, and I'm wondering what I can do to make it more efficient? What I've got so far is below. I wish to still record all steps. I think AppendTo is doing the most damage to the time, is there a faster alternative?
rk4[f_, variables_, valtinit_, tinit_, tfinal_, nsteps_] :=
Module[table, xlist, ylist, step, k1, k2, k3, k4,
xlist = tinit;
step = N[(tfinal - tinit)/(nsteps)];
ylist = valtinit;
table = xlist, ylist;
Table[
k1 = step* f /. MapThread[Rule, variables, ylist]; (*
Equivalent to step* f/.Thread[Rule[variables,ylist]]*)
k2 = step*f /. MapThread[Rule, variables, k1/2 + ylist];
k3 = step*f /. MapThread[Rule, variables, k2/2 + ylist];
k4 = step*f /. MapThread[Rule, variables, k3 + ylist];
ylist += 1/6 (k1 + 2 (k2 + k3) + k4);
xlist += step;
AppendTo[table, xlist, ylist];
xlist, ylist, nsteps];
table
];
Example Input:
funclist = -x + y, x - y;
initials = 1, 2;
variables = x, y;
init = 0;
final = 200;
nstep = 20000;
approx = rk4[funclist, variables, initials, init, final, nstep]//AbsoluteTiming;
3.59932,...
I'd love some suggestions!
differential-equations numerical-integration
asked 4 hours ago
ShinaolordShinaolord
808
$endgroup$
1
$begingroup$
AppendTois quadratic time complexity. Might be better to preallocate and set by index. Also it'll be much faster to not useRuleand instead code stuff up a little bit more explicitly. As a general rule, too, use vectorized operators. Those can be very fast. And if everything can be totally functional over "packed arrays" (look them up here) it'll be very quick too.
$endgroup$
– b3m2a1
3 hours ago
$begingroup$
I'll work on implementing it more explicity, this is what came to find first though. It'll require some changes to the inputs, I'll have to ponder this. And preallocating the list is a quick change that won't be an issue to do, I can't believe I forgot that's faster :(. Thanks though!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
Shinaoloard, usingJoin[ xlist, ylist, Table[ k1 = step*f /. MapThread[Rule, variables, ylist]; k2 = step*f /. MapThread[Rule, variables, k1/2 + ylist]; k3 = step*f /. MapThread[Rule, variables, k2/2 + ylist]; k4 = step*f /. MapThread[Rule, variables, k3 + ylist]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; xlist, ylist, nsteps ] ]as return value is already a first step. There is no point in appending if you use aTableanyways.
$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
@HenrikSchumacher do you think it would be faster to Pre-allocate a list of length nsteps, and append the values, or to join the values using table? I can obviously changeTabletoDoto remove the time it takes to make the table list, going by b3m2a1's method, or I could useJoinas you have suggested. I'm thinking your method may be faster, though. I've already removed theMapThreadpart, I am testing the speed increase granted by that at the moment. Just curious which path you think will be faster.
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
I am currently testing the speed difference between the one in the post andrk4t2[f_, valtinit_, tinit_, tfinal_, nsteps_] := Module[test, table, xlist, ylist, step, k1, k2, k3, k4, xlist = tinit; step = N[(tfinal - tinit)/(nsteps)]; ylist = valtinit; table = xlist, ylist; test = Table[ k1 = step* f[ylist] ; k2 = step*f[k1/2 + ylist]; k3 = step*f[k2/2 + ylist]; k4 = step*f[k3 + ylist]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; xlist, ylist, nsteps]; Join[table, test] ];
$endgroup$
– Shinaolord
3 hours ago
|
show 1 more comment
3
3
3
1
$begingroup$
So, I've implemented RK4, and I'm wondering what I can do to make it more efficient? What I've got so far is below. I wish to still record all steps. I think AppendTo is doing the most damage to the time, is there a faster alternative?
rk4[f_, variables_, valtinit_, tinit_, tfinal_, nsteps_] :=
Module[table, xlist, ylist, step, k1, k2, k3, k4,
xlist = tinit;
step = N[(tfinal - tinit)/(nsteps)];
ylist = valtinit;
table = xlist, ylist;
Table[
k1 = step* f /. MapThread[Rule, variables, ylist]; (*
Equivalent to step* f/.Thread[Rule[variables,ylist]]*)
k2 = step*f /. MapThread[Rule, variables, k1/2 + ylist];
k3 = step*f /. MapThread[Rule, variables, k2/2 + ylist];
k4 = step*f /. MapThread[Rule, variables, k3 + ylist];
ylist += 1/6 (k1 + 2 (k2 + k3) + k4);
xlist += step;
AppendTo[table, xlist, ylist];
xlist, ylist, nsteps];
table
];
Example Input:
funclist = -x + y, x - y;
initials = 1, 2;
variables = x, y;
init = 0;
final = 200;
nstep = 20000;
approx = rk4[funclist, variables, initials, init, final, nstep]//AbsoluteTiming;
3.59932,...
I'd love some suggestions!
differential-equations numerical-integration
asked 4 hours ago
ShinaolordShinaolord
808
$endgroup$
So, I've implemented RK4, and I'm wondering what I can do to make it more efficient? What I've got so far is below. I wish to still record all steps. I think AppendTo is doing the most damage to the time, is there a faster alternative?
rk4[f_, variables_, valtinit_, tinit_, tfinal_, nsteps_] :=
Module[table, xlist, ylist, step, k1, k2, k3, k4,
xlist = tinit;
step = N[(tfinal - tinit)/(nsteps)];
ylist = valtinit;
table = xlist, ylist;
Table[
k1 = step* f /. MapThread[Rule, variables, ylist]; (*
Equivalent to step* f/.Thread[Rule[variables,ylist]]*)
k2 = step*f /. MapThread[Rule, variables, k1/2 + ylist];
k3 = step*f /. MapThread[Rule, variables, k2/2 + ylist];
k4 = step*f /. MapThread[Rule, variables, k3 + ylist];
ylist += 1/6 (k1 + 2 (k2 + k3) + k4);
xlist += step;
AppendTo[table, xlist, ylist];
xlist, ylist, nsteps];
table
];
Example Input:
funclist = -x + y, x - y;
initials = 1, 2;
variables = x, y;
init = 0;
final = 200;
nstep = 20000;
approx = rk4[funclist, variables, initials, init, final, nstep]//AbsoluteTiming;
3.59932,...
I'd love some suggestions!
differential-equations numerical-integration
differential-equations numerical-integration
asked 4 hours ago
ShinaolordShinaolord
808
asked 4 hours ago
ShinaolordShinaolord
808
asked 4 hours ago
ShinaolordShinaolord
808
asked 4 hours ago
ShinaolordShinaolord
808
asked 4 hours ago
ShinaolordShinaolord
808
808
1
$begingroup$
AppendTois quadratic time complexity. Might be better to preallocate and set by index. Also it'll be much faster to not useRuleand instead code stuff up a little bit more explicitly. As a general rule, too, use vectorized operators. Those can be very fast. And if everything can be totally functional over "packed arrays" (look them up here) it'll be very quick too.
$endgroup$
– b3m2a1
3 hours ago
$begingroup$
I'll work on implementing it more explicity, this is what came to find first though. It'll require some changes to the inputs, I'll have to ponder this. And preallocating the list is a quick change that won't be an issue to do, I can't believe I forgot that's faster :(. Thanks though!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
Shinaoloard, usingJoin[ xlist, ylist, Table[ k1 = step*f /. MapThread[Rule, variables, ylist]; k2 = step*f /. MapThread[Rule, variables, k1/2 + ylist]; k3 = step*f /. MapThread[Rule, variables, k2/2 + ylist]; k4 = step*f /. MapThread[Rule, variables, k3 + ylist]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; xlist, ylist, nsteps ] ]as return value is already a first step. There is no point in appending if you use aTableanyways.
$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
@HenrikSchumacher do you think it would be faster to Pre-allocate a list of length nsteps, and append the values, or to join the values using table? I can obviously changeTabletoDoto remove the time it takes to make the table list, going by b3m2a1's method, or I could useJoinas you have suggested. I'm thinking your method may be faster, though. I've already removed theMapThreadpart, I am testing the speed increase granted by that at the moment. Just curious which path you think will be faster.
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
I am currently testing the speed difference between the one in the post andrk4t2[f_, valtinit_, tinit_, tfinal_, nsteps_] := Module[test, table, xlist, ylist, step, k1, k2, k3, k4, xlist = tinit; step = N[(tfinal - tinit)/(nsteps)]; ylist = valtinit; table = xlist, ylist; test = Table[ k1 = step* f[ylist] ; k2 = step*f[k1/2 + ylist]; k3 = step*f[k2/2 + ylist]; k4 = step*f[k3 + ylist]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; xlist, ylist, nsteps]; Join[table, test] ];
$endgroup$
– Shinaolord
3 hours ago
|
show 1 more comment
1
$begingroup$
AppendTois quadratic time complexity. Might be better to preallocate and set by index. Also it'll be much faster to not useRuleand instead code stuff up a little bit more explicitly. As a general rule, too, use vectorized operators. Those can be very fast. And if everything can be totally functional over "packed arrays" (look them up here) it'll be very quick too.
$endgroup$
– b3m2a1
3 hours ago
$begingroup$
I'll work on implementing it more explicity, this is what came to find first though. It'll require some changes to the inputs, I'll have to ponder this. And preallocating the list is a quick change that won't be an issue to do, I can't believe I forgot that's faster :(. Thanks though!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
Shinaoloard, usingJoin[ xlist, ylist, Table[ k1 = step*f /. MapThread[Rule, variables, ylist]; k2 = step*f /. MapThread[Rule, variables, k1/2 + ylist]; k3 = step*f /. MapThread[Rule, variables, k2/2 + ylist]; k4 = step*f /. MapThread[Rule, variables, k3 + ylist]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; xlist, ylist, nsteps ] ]as return value is already a first step. There is no point in appending if you use aTableanyways.
$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
@HenrikSchumacher do you think it would be faster to Pre-allocate a list of length nsteps, and append the values, or to join the values using table? I can obviously changeTabletoDoto remove the time it takes to make the table list, going by b3m2a1's method, or I could useJoinas you have suggested. I'm thinking your method may be faster, though. I've already removed theMapThreadpart, I am testing the speed increase granted by that at the moment. Just curious which path you think will be faster.
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
I am currently testing the speed difference between the one in the post andrk4t2[f_, valtinit_, tinit_, tfinal_, nsteps_] := Module[test, table, xlist, ylist, step, k1, k2, k3, k4, xlist = tinit; step = N[(tfinal - tinit)/(nsteps)]; ylist = valtinit; table = xlist, ylist; test = Table[ k1 = step* f[ylist] ; k2 = step*f[k1/2 + ylist]; k3 = step*f[k2/2 + ylist]; k4 = step*f[k3 + ylist]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; xlist, ylist, nsteps]; Join[table, test] ];
$endgroup$
– Shinaolord
3 hours ago
1
1
$begingroup$
AppendTo is quadratic time complexity. Might be better to preallocate and set by index. Also it'll be much faster to not use Rule and instead code stuff up a little bit more explicitly. As a general rule, too, use vectorized operators. Those can be very fast. And if everything can be totally functional over "packed arrays" (look them up here) it'll be very quick too.$endgroup$
– b3m2a1
3 hours ago
$begingroup$
AppendTo is quadratic time complexity. Might be better to preallocate and set by index. Also it'll be much faster to not use Rule and instead code stuff up a little bit more explicitly. As a general rule, too, use vectorized operators. Those can be very fast. And if everything can be totally functional over "packed arrays" (look them up here) it'll be very quick too.$endgroup$
– b3m2a1
3 hours ago
$begingroup$
I'll work on implementing it more explicity, this is what came to find first though. It'll require some changes to the inputs, I'll have to ponder this. And preallocating the list is a quick change that won't be an issue to do, I can't believe I forgot that's faster :(. Thanks though!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
I'll work on implementing it more explicity, this is what came to find first though. It'll require some changes to the inputs, I'll have to ponder this. And preallocating the list is a quick change that won't be an issue to do, I can't believe I forgot that's faster :(. Thanks though!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
Shinaoloard, using
Join[ xlist, ylist, Table[ k1 = step*f /. MapThread[Rule, variables, ylist]; k2 = step*f /. MapThread[Rule, variables, k1/2 + ylist]; k3 = step*f /. MapThread[Rule, variables, k2/2 + ylist]; k4 = step*f /. MapThread[Rule, variables, k3 + ylist]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; xlist, ylist, nsteps ] ] as return value is already a first step. There is no point in appending if you use a Table anyways.$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
Shinaoloard, using
Join[ xlist, ylist, Table[ k1 = step*f /. MapThread[Rule, variables, ylist]; k2 = step*f /. MapThread[Rule, variables, k1/2 + ylist]; k3 = step*f /. MapThread[Rule, variables, k2/2 + ylist]; k4 = step*f /. MapThread[Rule, variables, k3 + ylist]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; xlist, ylist, nsteps ] ] as return value is already a first step. There is no point in appending if you use a Table anyways.$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
@HenrikSchumacher do you think it would be faster to Pre-allocate a list of length nsteps, and append the values, or to join the values using table? I can obviously change
Table to Do to remove the time it takes to make the table list, going by b3m2a1's method, or I could use Join as you have suggested. I'm thinking your method may be faster, though. I've already removed the MapThread part, I am testing the speed increase granted by that at the moment. Just curious which path you think will be faster.$endgroup$
– Shinaolord
3 hours ago
$begingroup$
@HenrikSchumacher do you think it would be faster to Pre-allocate a list of length nsteps, and append the values, or to join the values using table? I can obviously change
Table to Do to remove the time it takes to make the table list, going by b3m2a1's method, or I could use Join as you have suggested. I'm thinking your method may be faster, though. I've already removed the MapThread part, I am testing the speed increase granted by that at the moment. Just curious which path you think will be faster.$endgroup$
– Shinaolord
3 hours ago
$begingroup$
I am currently testing the speed difference between the one in the post and
rk4t2[f_, valtinit_, tinit_, tfinal_, nsteps_] := Module[test, table, xlist, ylist, step, k1, k2, k3, k4, xlist = tinit; step = N[(tfinal - tinit)/(nsteps)]; ylist = valtinit; table = xlist, ylist; test = Table[ k1 = step* f[ylist] ; k2 = step*f[k1/2 + ylist]; k3 = step*f[k2/2 + ylist]; k4 = step*f[k3 + ylist]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; xlist, ylist, nsteps]; Join[table, test] ];$endgroup$
– Shinaolord
3 hours ago
$begingroup$
I am currently testing the speed difference between the one in the post and
rk4t2[f_, valtinit_, tinit_, tfinal_, nsteps_] := Module[test, table, xlist, ylist, step, k1, k2, k3, k4, xlist = tinit; step = N[(tfinal - tinit)/(nsteps)]; ylist = valtinit; table = xlist, ylist; test = Table[ k1 = step* f[ylist] ; k2 = step*f[k1/2 + ylist]; k3 = step*f[k2/2 + ylist]; k4 = step*f[k3 + ylist]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; xlist, ylist, nsteps]; Join[table, test] ];$endgroup$
– Shinaolord
3 hours ago
|
show 1 more comment
1 Answer
1
active
oldest
votes
5
$begingroup$
Just to give you an impression how fast things may get when you use the right tools.
For given stepsize τ and given vectorfield F, this creates a CompiledFunction cStep that computes a single Runge-Kutta step
F = X [Function] -Indexed[X, 2], Indexed[X, 1];
τ = 0.01;
Block[YY, Y, k1, k2, k3, k4,
YY = Table[Compile`GetElement[Y, i], i, 1, 2];
k1 = τ F[YY];
k2 = τ F[0.5 k1 + YY];
k3 = τ F[0.5 k2 + YY];
k4 = τ F[k3 + YY];
cStep = With[code = YY + (k1 + 2. (k2 + k3) + k4)/6. ,
Compile[Y, _Real, 1,
code,
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
]
]
];
Now we can apply it 2 million times with NestList and still need only 2 seconds.
nsteps = 20000000;
xlist = Range[0., step nsteps, step];
Ylist = NestList[cStep, initials, nsteps]; // AbsoluteTiming // First
2.08678
answered 3 hours ago
Henrik SchumacherHenrik Schumacher
57.9k579159
$endgroup$
$begingroup$
Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
You're welcome.
$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
I'll have to play around withCompile, it definitely seems like a massive speed up if used correctly.
$endgroup$
– Shinaolord
3 hours ago
add a 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: "387"
;
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: false,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
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
);
);
draft saved
draft discarded
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
var $window = $(window),
onScroll = function(e)
var $elem = $('.new-login-left'),
docViewTop = $window.scrollTop(),
docViewBottom = docViewTop + $window.height(),
elemTop = $elem.offset().top,
elemBottom = elemTop + $elem.height();
if ((docViewTop elemBottom))
StackExchange.using('gps', function() StackExchange.gps.track('embedded_signup_form.view', location: 'question_page' ); );
$window.unbind('scroll', onScroll);
;
$window.on('scroll', onScroll);
);
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%2fmathematica.stackexchange.com%2fquestions%2f194002%2fways-to-speed-up-user-implemented-rk4%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
5
$begingroup$
Just to give you an impression how fast things may get when you use the right tools.
For given stepsize τ and given vectorfield F, this creates a CompiledFunction cStep that computes a single Runge-Kutta step
F = X [Function] -Indexed[X, 2], Indexed[X, 1];
τ = 0.01;
Block[YY, Y, k1, k2, k3, k4,
YY = Table[Compile`GetElement[Y, i], i, 1, 2];
k1 = τ F[YY];
k2 = τ F[0.5 k1 + YY];
k3 = τ F[0.5 k2 + YY];
k4 = τ F[k3 + YY];
cStep = With[code = YY + (k1 + 2. (k2 + k3) + k4)/6. ,
Compile[Y, _Real, 1,
code,
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
]
]
];
Now we can apply it 2 million times with NestList and still need only 2 seconds.
nsteps = 20000000;
xlist = Range[0., step nsteps, step];
Ylist = NestList[cStep, initials, nsteps]; // AbsoluteTiming // First
2.08678
answered 3 hours ago
Henrik SchumacherHenrik Schumacher
57.9k579159
$endgroup$
$begingroup$
Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
You're welcome.
$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
I'll have to play around withCompile, it definitely seems like a massive speed up if used correctly.
$endgroup$
– Shinaolord
3 hours ago
add a comment |
5
$begingroup$
Just to give you an impression how fast things may get when you use the right tools.
For given stepsize τ and given vectorfield F, this creates a CompiledFunction cStep that computes a single Runge-Kutta step
F = X [Function] -Indexed[X, 2], Indexed[X, 1];
τ = 0.01;
Block[YY, Y, k1, k2, k3, k4,
YY = Table[Compile`GetElement[Y, i], i, 1, 2];
k1 = τ F[YY];
k2 = τ F[0.5 k1 + YY];
k3 = τ F[0.5 k2 + YY];
k4 = τ F[k3 + YY];
cStep = With[code = YY + (k1 + 2. (k2 + k3) + k4)/6. ,
Compile[Y, _Real, 1,
code,
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
]
]
];
Now we can apply it 2 million times with NestList and still need only 2 seconds.
nsteps = 20000000;
xlist = Range[0., step nsteps, step];
Ylist = NestList[cStep, initials, nsteps]; // AbsoluteTiming // First
2.08678
answered 3 hours ago
Henrik SchumacherHenrik Schumacher
57.9k579159
$endgroup$
$begingroup$
Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
You're welcome.
$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
I'll have to play around withCompile, it definitely seems like a massive speed up if used correctly.
$endgroup$
– Shinaolord
3 hours ago
add a comment |
5
5
5
$begingroup$
Just to give you an impression how fast things may get when you use the right tools.
For given stepsize τ and given vectorfield F, this creates a CompiledFunction cStep that computes a single Runge-Kutta step
F = X [Function] -Indexed[X, 2], Indexed[X, 1];
τ = 0.01;
Block[YY, Y, k1, k2, k3, k4,
YY = Table[Compile`GetElement[Y, i], i, 1, 2];
k1 = τ F[YY];
k2 = τ F[0.5 k1 + YY];
k3 = τ F[0.5 k2 + YY];
k4 = τ F[k3 + YY];
cStep = With[code = YY + (k1 + 2. (k2 + k3) + k4)/6. ,
Compile[Y, _Real, 1,
code,
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
]
]
];
Now we can apply it 2 million times with NestList and still need only 2 seconds.
nsteps = 20000000;
xlist = Range[0., step nsteps, step];
Ylist = NestList[cStep, initials, nsteps]; // AbsoluteTiming // First
2.08678
answered 3 hours ago
Henrik SchumacherHenrik Schumacher
57.9k579159
$endgroup$
Just to give you an impression how fast things may get when you use the right tools.
For given stepsize τ and given vectorfield F, this creates a CompiledFunction cStep that computes a single Runge-Kutta step
F = X [Function] -Indexed[X, 2], Indexed[X, 1];
τ = 0.01;
Block[YY, Y, k1, k2, k3, k4,
YY = Table[Compile`GetElement[Y, i], i, 1, 2];
k1 = τ F[YY];
k2 = τ F[0.5 k1 + YY];
k3 = τ F[0.5 k2 + YY];
k4 = τ F[k3 + YY];
cStep = With[code = YY + (k1 + 2. (k2 + k3) + k4)/6. ,
Compile[Y, _Real, 1,
code,
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
]
]
];
Now we can apply it 2 million times with NestList and still need only 2 seconds.
nsteps = 20000000;
xlist = Range[0., step nsteps, step];
Ylist = NestList[cStep, initials, nsteps]; // AbsoluteTiming // First
2.08678
answered 3 hours ago
Henrik SchumacherHenrik Schumacher
57.9k579159
answered 3 hours ago
Henrik SchumacherHenrik Schumacher
57.9k579159
answered 3 hours ago
Henrik SchumacherHenrik Schumacher
57.9k579159
answered 3 hours ago
Henrik SchumacherHenrik Schumacher
57.9k579159
57.9k579159
$begingroup$
Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
You're welcome.
$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
I'll have to play around withCompile, it definitely seems like a massive speed up if used correctly.
$endgroup$
– Shinaolord
3 hours ago
add a comment |
$begingroup$
Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
You're welcome.
$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
I'll have to play around withCompile, it definitely seems like a massive speed up if used correctly.
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
You're welcome.
$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
You're welcome.
$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
I'll have to play around with
Compile, it definitely seems like a massive speed up if used correctly.$endgroup$
– Shinaolord
3 hours ago
$begingroup$
I'll have to play around with
Compile, it definitely seems like a massive speed up if used correctly.$endgroup$
– Shinaolord
3 hours ago
add a comment |
draft saved
draft discarded
Thanks for contributing an answer to Mathematica 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');
var $window = $(window),
onScroll = function(e)
var $elem = $('.new-login-left'),
docViewTop = $window.scrollTop(),
docViewBottom = docViewTop + $window.height(),
elemTop = $elem.offset().top,
elemBottom = elemTop + $elem.height();
if ((docViewTop elemBottom))
StackExchange.using('gps', function() StackExchange.gps.track('embedded_signup_form.view', location: 'question_page' ); );
$window.unbind('scroll', onScroll);
;
$window.on('scroll', onScroll);
);
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%2fmathematica.stackexchange.com%2fquestions%2f194002%2fways-to-speed-up-user-implemented-rk4%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');
var $window = $(window),
onScroll = function(e)
var $elem = $('.new-login-left'),
docViewTop = $window.scrollTop(),
docViewBottom = docViewTop + $window.height(),
elemTop = $elem.offset().top,
elemBottom = elemTop + $elem.height();
if ((docViewTop elemBottom))
StackExchange.using('gps', function() StackExchange.gps.track('embedded_signup_form.view', location: 'question_page' ); );
$window.unbind('scroll', onScroll);
;
$window.on('scroll', onScroll);
);
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');
var $window = $(window),
onScroll = function(e)
var $elem = $('.new-login-left'),
docViewTop = $window.scrollTop(),
docViewBottom = docViewTop + $window.height(),
elemTop = $elem.offset().top,
elemBottom = elemTop + $elem.height();
if ((docViewTop elemBottom))
StackExchange.using('gps', function() StackExchange.gps.track('embedded_signup_form.view', location: 'question_page' ); );
$window.unbind('scroll', onScroll);
;
$window.on('scroll', onScroll);
);
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');
var $window = $(window),
onScroll = function(e)
var $elem = $('.new-login-left'),
docViewTop = $window.scrollTop(),
docViewBottom = docViewTop + $window.height(),
elemTop = $elem.offset().top,
elemBottom = elemTop + $elem.height();
if ((docViewTop elemBottom))
StackExchange.using('gps', function() StackExchange.gps.track('embedded_signup_form.view', location: 'question_page' ); );
$window.unbind('scroll', onScroll);
;
$window.on('scroll', onScroll);
);
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
1
$begingroup$
AppendTois quadratic time complexity. Might be better to preallocate and set by index. Also it'll be much faster to not useRuleand instead code stuff up a little bit more explicitly. As a general rule, too, use vectorized operators. Those can be very fast. And if everything can be totally functional over "packed arrays" (look them up here) it'll be very quick too.$endgroup$
– b3m2a1
3 hours ago
$begingroup$
I'll work on implementing it more explicity, this is what came to find first though. It'll require some changes to the inputs, I'll have to ponder this. And preallocating the list is a quick change that won't be an issue to do, I can't believe I forgot that's faster :(. Thanks though!
$endgroup$
– Shinaolord
3 hours ago
$begingroup$
Shinaoloard, using
Join[ xlist, ylist, Table[ k1 = step*f /. MapThread[Rule, variables, ylist]; k2 = step*f /. MapThread[Rule, variables, k1/2 + ylist]; k3 = step*f /. MapThread[Rule, variables, k2/2 + ylist]; k4 = step*f /. MapThread[Rule, variables, k3 + ylist]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; xlist, ylist, nsteps ] ]as return value is already a first step. There is no point in appending if you use aTableanyways.$endgroup$
– Henrik Schumacher
3 hours ago
$begingroup$
@HenrikSchumacher do you think it would be faster to Pre-allocate a list of length nsteps, and append the values, or to join the values using table? I can obviously change
TabletoDoto remove the time it takes to make the table list, going by b3m2a1's method, or I could useJoinas you have suggested. I'm thinking your method may be faster, though. I've already removed theMapThreadpart, I am testing the speed increase granted by that at the moment. Just curious which path you think will be faster.$endgroup$
– Shinaolord
3 hours ago
$begingroup$
I am currently testing the speed difference between the one in the post and
rk4t2[f_, valtinit_, tinit_, tfinal_, nsteps_] := Module[test, table, xlist, ylist, step, k1, k2, k3, k4, xlist = tinit; step = N[(tfinal - tinit)/(nsteps)]; ylist = valtinit; table = xlist, ylist; test = Table[ k1 = step* f[ylist] ; k2 = step*f[k1/2 + ylist]; k3 = step*f[k2/2 + ylist]; k4 = step*f[k3 + ylist]; ylist += 1/6 (k1 + 2 (k2 + k3) + k4); xlist += step; xlist, ylist, nsteps]; Join[table, test] ];$endgroup$
– Shinaolord
3 hours ago