Ways to speed up user implemented RK4 The Next CEO of Stack OverflowSpeed up Numerical IntegrationHow to choose MaxStepFraction for optimal speed of NDSolveNIntegrate: how to speed up code?Compiling FoldList implementation for RK4How to speed up this code?Solving an unstable BVP numerically, accurately and efficientlyHow to speed up integral of results of PDE modelSolve BVP involving user defined functionWays to speed up PickSpeed up ParametricNDSolve

Can you teleport closer to a creature you are Frightened of?

How to compactly explain secondary and tertiary characters without resorting to stereotypes?

Early programmable calculators with RS-232

What difference does it make matching a word with/without a trailing whitespace?

Does the Idaho Potato Commission associate potato skins with healthy eating?

Ising model simulation

Masking layers by a vector polygon layer in QGIS

Do I need to write [sic] when including a quotation with a number less than 10 that isn't written out?

Is a linearly independent set whose span is dense a Schauder basis?

"Eavesdropping" vs "Listen in on"

Does int main() need a declaration on C++?

How can I replace x-axis labels with pre-determined symbols?

Could you use a laser beam as a modulated carrier wave for radio signal?

Oldie but Goldie

Compensation for working overtime on Saturdays

Upgrading From a 9 Speed Sora Derailleur?

How to show a landlord what we have in savings?

My ex-girlfriend uses my Apple ID to login to her iPad, do I have to give her my Apple ID password to reset it?

Find a path from s to t using as few red nodes as possible

Is this a new Fibonacci Identity?

What did the word "leisure" mean in late 18th Century usage?

Physiological effects of huge anime eyes

Planeswalker Ability and Death Timing

How to coordinate airplane tickets?



Ways to speed up user implemented RK4



The Next CEO of Stack OverflowSpeed up Numerical IntegrationHow to choose MaxStepFraction for optimal speed of NDSolveNIntegrate: how to speed up code?Compiling FoldList implementation for RK4How to speed up this code?Solving an unstable BVP numerically, accurately and efficientlyHow to speed up integral of results of PDE modelSolve BVP involving user defined functionWays to speed up PickSpeed up ParametricNDSolve










8












$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!










share|improve this question











$endgroup$







  • 3




    $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
    Mar 26 at 21:31










  • $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
    Mar 26 at 21:33










  • $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
    Mar 26 at 21:40






  • 3




    $begingroup$
    Why not just get NDSolve[] to use fourth-order Runge-Kutta to begin with?
    $endgroup$
    – J. M. is slightly pensive
    Mar 27 at 2:46






  • 1




    $begingroup$
    @J.M.isslightlypensive I know it can, I just wanted to make sure I could actually code it myself, instead of just using options to get Mathematica to do it for me (:. Thanks for trying to help though!!
    $endgroup$
    – Shinaolord
    Mar 27 at 13:27















8












$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!










share|improve this question











$endgroup$







  • 3




    $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
    Mar 26 at 21:31










  • $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
    Mar 26 at 21:33










  • $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
    Mar 26 at 21:40






  • 3




    $begingroup$
    Why not just get NDSolve[] to use fourth-order Runge-Kutta to begin with?
    $endgroup$
    – J. M. is slightly pensive
    Mar 27 at 2:46






  • 1




    $begingroup$
    @J.M.isslightlypensive I know it can, I just wanted to make sure I could actually code it myself, instead of just using options to get Mathematica to do it for me (:. Thanks for trying to help though!!
    $endgroup$
    – Shinaolord
    Mar 27 at 13:27













8












8








8


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!










share|improve this question











$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 performance-tuning






share|improve this question















share|improve this question













share|improve this question




share|improve this question








edited Mar 27 at 2:17









xzczd

27.4k573255




27.4k573255










asked Mar 26 at 21:21









ShinaolordShinaolord

1088




1088







  • 3




    $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
    Mar 26 at 21:31










  • $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
    Mar 26 at 21:33










  • $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
    Mar 26 at 21:40






  • 3




    $begingroup$
    Why not just get NDSolve[] to use fourth-order Runge-Kutta to begin with?
    $endgroup$
    – J. M. is slightly pensive
    Mar 27 at 2:46






  • 1




    $begingroup$
    @J.M.isslightlypensive I know it can, I just wanted to make sure I could actually code it myself, instead of just using options to get Mathematica to do it for me (:. Thanks for trying to help though!!
    $endgroup$
    – Shinaolord
    Mar 27 at 13:27












  • 3




    $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
    Mar 26 at 21:31










  • $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
    Mar 26 at 21:33










  • $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
    Mar 26 at 21:40






  • 3




    $begingroup$
    Why not just get NDSolve[] to use fourth-order Runge-Kutta to begin with?
    $endgroup$
    – J. M. is slightly pensive
    Mar 27 at 2:46






  • 1




    $begingroup$
    @J.M.isslightlypensive I know it can, I just wanted to make sure I could actually code it myself, instead of just using options to get Mathematica to do it for me (:. Thanks for trying to help though!!
    $endgroup$
    – Shinaolord
    Mar 27 at 13:27







3




3




$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
Mar 26 at 21:31




$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
Mar 26 at 21:31












$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
Mar 26 at 21:33




$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
Mar 26 at 21:33












$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
Mar 26 at 21:40




$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
Mar 26 at 21:40




3




3




$begingroup$
Why not just get NDSolve[] to use fourth-order Runge-Kutta to begin with?
$endgroup$
– J. M. is slightly pensive
Mar 27 at 2:46




$begingroup$
Why not just get NDSolve[] to use fourth-order Runge-Kutta to begin with?
$endgroup$
– J. M. is slightly pensive
Mar 27 at 2:46




1




1




$begingroup$
@J.M.isslightlypensive I know it can, I just wanted to make sure I could actually code it myself, instead of just using options to get Mathematica to do it for me (:. Thanks for trying to help though!!
$endgroup$
– Shinaolord
Mar 27 at 13:27




$begingroup$
@J.M.isslightlypensive I know it can, I just wanted to make sure I could actually code it myself, instead of just using options to get Mathematica to do it for me (:. Thanks for trying to help though!!
$endgroup$
– Shinaolord
Mar 27 at 13:27










1 Answer
1






active

oldest

votes


















16












$begingroup$

Just to give you an impression how fast things may get when you use the right tools.



For given stepsize τ and given vector field 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 20 million times with NestList and it stills takes only about 2 seconds.



nsteps = 20000000;
xlist = Range[0., τ nsteps, τ];
Ylist = NestList[cStep, 1., 0., nsteps]; // AbsoluteTiming // First



2.08678




Edit



This can be sped up even more my avoiding NestList (the loop behind it can also be compiled which saves several calls to libraries) and by utilizing that the dimension of the ODE is known at compile time. For low dimensional systems, it may be also beneficial to avoid tensor operations altogether and to perform computations in scalar registers as done below.



τ = 0.01;
cFlow = Block[YY, Y, k1, k2, k3, k4, τ, Ylist, j,
YY = Table[Compile`GetElement[Ylist, j, i], i, 1, 2];
k1 = τ F[YY];
k2 = τ F[0.5 k1 + YY];
k3 = τ F[0.5 k2 + YY];
k4 = τ F[k3 + YY];
With[
code1 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[1]],
code2 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[2]]
,
Compile[Y0, _Real, 1, τ, _Real, n, _Integer,
Block[Ylist,
Ylist = Table[0., n + 1, Length[Y0]];
Ylist[[1]] = Y0;
Do[
Ylist[[j + 1, 1]] = code1;
Ylist[[j + 1, 2]] = code2;
,
j, 1, n];
Ylist
],
CompilationTarget -> "C", RuntimeOptions -> "Speed"
]
]
];
Ylist2 = cFlow[1., 0., τ, nsteps]; // AbsoluteTiming // First



1.06549




Don't be too upset by parts of the code being highlighted in red; this is on purpose.






share|improve this answer











$endgroup$








  • 1




    $begingroup$
    Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
    $endgroup$
    – Shinaolord
    Mar 26 at 22:07










  • $begingroup$
    You're welcome.
    $endgroup$
    – Henrik Schumacher
    Mar 26 at 22:08










  • $begingroup$
    I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
    $endgroup$
    – Shinaolord
    Mar 26 at 22:08






  • 1




    $begingroup$
    This is exactly the kind of thing I like to show people when they complain about the slowness of Mathematica. Of course with some cleverness in vectorized operations Compile could probably be avoided altogether if one so desired.
    $endgroup$
    – b3m2a1
    Mar 27 at 7:45











  • $begingroup$
    @b3m2a1 Yeah, right. However, with ODE systems of this small dimension (dim = 2) I am not sure how to utilize vectorization since the ODE has to be solved in serial and because the vector field F may be very nonlinear.
    $endgroup$
    – Henrik Schumacher
    Mar 27 at 7:53











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


















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









16












$begingroup$

Just to give you an impression how fast things may get when you use the right tools.



For given stepsize τ and given vector field 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 20 million times with NestList and it stills takes only about 2 seconds.



nsteps = 20000000;
xlist = Range[0., τ nsteps, τ];
Ylist = NestList[cStep, 1., 0., nsteps]; // AbsoluteTiming // First



2.08678




Edit



This can be sped up even more my avoiding NestList (the loop behind it can also be compiled which saves several calls to libraries) and by utilizing that the dimension of the ODE is known at compile time. For low dimensional systems, it may be also beneficial to avoid tensor operations altogether and to perform computations in scalar registers as done below.



τ = 0.01;
cFlow = Block[YY, Y, k1, k2, k3, k4, τ, Ylist, j,
YY = Table[Compile`GetElement[Ylist, j, i], i, 1, 2];
k1 = τ F[YY];
k2 = τ F[0.5 k1 + YY];
k3 = τ F[0.5 k2 + YY];
k4 = τ F[k3 + YY];
With[
code1 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[1]],
code2 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[2]]
,
Compile[Y0, _Real, 1, τ, _Real, n, _Integer,
Block[Ylist,
Ylist = Table[0., n + 1, Length[Y0]];
Ylist[[1]] = Y0;
Do[
Ylist[[j + 1, 1]] = code1;
Ylist[[j + 1, 2]] = code2;
,
j, 1, n];
Ylist
],
CompilationTarget -> "C", RuntimeOptions -> "Speed"
]
]
];
Ylist2 = cFlow[1., 0., τ, nsteps]; // AbsoluteTiming // First



1.06549




Don't be too upset by parts of the code being highlighted in red; this is on purpose.






share|improve this answer











$endgroup$








  • 1




    $begingroup$
    Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
    $endgroup$
    – Shinaolord
    Mar 26 at 22:07










  • $begingroup$
    You're welcome.
    $endgroup$
    – Henrik Schumacher
    Mar 26 at 22:08










  • $begingroup$
    I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
    $endgroup$
    – Shinaolord
    Mar 26 at 22:08






  • 1




    $begingroup$
    This is exactly the kind of thing I like to show people when they complain about the slowness of Mathematica. Of course with some cleverness in vectorized operations Compile could probably be avoided altogether if one so desired.
    $endgroup$
    – b3m2a1
    Mar 27 at 7:45











  • $begingroup$
    @b3m2a1 Yeah, right. However, with ODE systems of this small dimension (dim = 2) I am not sure how to utilize vectorization since the ODE has to be solved in serial and because the vector field F may be very nonlinear.
    $endgroup$
    – Henrik Schumacher
    Mar 27 at 7:53















16












$begingroup$

Just to give you an impression how fast things may get when you use the right tools.



For given stepsize τ and given vector field 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 20 million times with NestList and it stills takes only about 2 seconds.



nsteps = 20000000;
xlist = Range[0., τ nsteps, τ];
Ylist = NestList[cStep, 1., 0., nsteps]; // AbsoluteTiming // First



2.08678




Edit



This can be sped up even more my avoiding NestList (the loop behind it can also be compiled which saves several calls to libraries) and by utilizing that the dimension of the ODE is known at compile time. For low dimensional systems, it may be also beneficial to avoid tensor operations altogether and to perform computations in scalar registers as done below.



τ = 0.01;
cFlow = Block[YY, Y, k1, k2, k3, k4, τ, Ylist, j,
YY = Table[Compile`GetElement[Ylist, j, i], i, 1, 2];
k1 = τ F[YY];
k2 = τ F[0.5 k1 + YY];
k3 = τ F[0.5 k2 + YY];
k4 = τ F[k3 + YY];
With[
code1 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[1]],
code2 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[2]]
,
Compile[Y0, _Real, 1, τ, _Real, n, _Integer,
Block[Ylist,
Ylist = Table[0., n + 1, Length[Y0]];
Ylist[[1]] = Y0;
Do[
Ylist[[j + 1, 1]] = code1;
Ylist[[j + 1, 2]] = code2;
,
j, 1, n];
Ylist
],
CompilationTarget -> "C", RuntimeOptions -> "Speed"
]
]
];
Ylist2 = cFlow[1., 0., τ, nsteps]; // AbsoluteTiming // First



1.06549




Don't be too upset by parts of the code being highlighted in red; this is on purpose.






share|improve this answer











$endgroup$








  • 1




    $begingroup$
    Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
    $endgroup$
    – Shinaolord
    Mar 26 at 22:07










  • $begingroup$
    You're welcome.
    $endgroup$
    – Henrik Schumacher
    Mar 26 at 22:08










  • $begingroup$
    I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
    $endgroup$
    – Shinaolord
    Mar 26 at 22:08






  • 1




    $begingroup$
    This is exactly the kind of thing I like to show people when they complain about the slowness of Mathematica. Of course with some cleverness in vectorized operations Compile could probably be avoided altogether if one so desired.
    $endgroup$
    – b3m2a1
    Mar 27 at 7:45











  • $begingroup$
    @b3m2a1 Yeah, right. However, with ODE systems of this small dimension (dim = 2) I am not sure how to utilize vectorization since the ODE has to be solved in serial and because the vector field F may be very nonlinear.
    $endgroup$
    – Henrik Schumacher
    Mar 27 at 7:53













16












16








16





$begingroup$

Just to give you an impression how fast things may get when you use the right tools.



For given stepsize τ and given vector field 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 20 million times with NestList and it stills takes only about 2 seconds.



nsteps = 20000000;
xlist = Range[0., τ nsteps, τ];
Ylist = NestList[cStep, 1., 0., nsteps]; // AbsoluteTiming // First



2.08678




Edit



This can be sped up even more my avoiding NestList (the loop behind it can also be compiled which saves several calls to libraries) and by utilizing that the dimension of the ODE is known at compile time. For low dimensional systems, it may be also beneficial to avoid tensor operations altogether and to perform computations in scalar registers as done below.



τ = 0.01;
cFlow = Block[YY, Y, k1, k2, k3, k4, τ, Ylist, j,
YY = Table[Compile`GetElement[Ylist, j, i], i, 1, 2];
k1 = τ F[YY];
k2 = τ F[0.5 k1 + YY];
k3 = τ F[0.5 k2 + YY];
k4 = τ F[k3 + YY];
With[
code1 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[1]],
code2 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[2]]
,
Compile[Y0, _Real, 1, τ, _Real, n, _Integer,
Block[Ylist,
Ylist = Table[0., n + 1, Length[Y0]];
Ylist[[1]] = Y0;
Do[
Ylist[[j + 1, 1]] = code1;
Ylist[[j + 1, 2]] = code2;
,
j, 1, n];
Ylist
],
CompilationTarget -> "C", RuntimeOptions -> "Speed"
]
]
];
Ylist2 = cFlow[1., 0., τ, nsteps]; // AbsoluteTiming // First



1.06549




Don't be too upset by parts of the code being highlighted in red; this is on purpose.






share|improve this answer











$endgroup$



Just to give you an impression how fast things may get when you use the right tools.



For given stepsize τ and given vector field 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 20 million times with NestList and it stills takes only about 2 seconds.



nsteps = 20000000;
xlist = Range[0., τ nsteps, τ];
Ylist = NestList[cStep, 1., 0., nsteps]; // AbsoluteTiming // First



2.08678




Edit



This can be sped up even more my avoiding NestList (the loop behind it can also be compiled which saves several calls to libraries) and by utilizing that the dimension of the ODE is known at compile time. For low dimensional systems, it may be also beneficial to avoid tensor operations altogether and to perform computations in scalar registers as done below.



τ = 0.01;
cFlow = Block[YY, Y, k1, k2, k3, k4, τ, Ylist, j,
YY = Table[Compile`GetElement[Ylist, j, i], i, 1, 2];
k1 = τ F[YY];
k2 = τ F[0.5 k1 + YY];
k3 = τ F[0.5 k2 + YY];
k4 = τ F[k3 + YY];
With[
code1 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[1]],
code2 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[2]]
,
Compile[Y0, _Real, 1, τ, _Real, n, _Integer,
Block[Ylist,
Ylist = Table[0., n + 1, Length[Y0]];
Ylist[[1]] = Y0;
Do[
Ylist[[j + 1, 1]] = code1;
Ylist[[j + 1, 2]] = code2;
,
j, 1, n];
Ylist
],
CompilationTarget -> "C", RuntimeOptions -> "Speed"
]
]
];
Ylist2 = cFlow[1., 0., τ, nsteps]; // AbsoluteTiming // First



1.06549




Don't be too upset by parts of the code being highlighted in red; this is on purpose.







share|improve this answer














share|improve this answer



share|improve this answer








edited Mar 27 at 14:19

























answered Mar 26 at 22:05









Henrik SchumacherHenrik Schumacher

58.7k581162




58.7k581162







  • 1




    $begingroup$
    Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
    $endgroup$
    – Shinaolord
    Mar 26 at 22:07










  • $begingroup$
    You're welcome.
    $endgroup$
    – Henrik Schumacher
    Mar 26 at 22:08










  • $begingroup$
    I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
    $endgroup$
    – Shinaolord
    Mar 26 at 22:08






  • 1




    $begingroup$
    This is exactly the kind of thing I like to show people when they complain about the slowness of Mathematica. Of course with some cleverness in vectorized operations Compile could probably be avoided altogether if one so desired.
    $endgroup$
    – b3m2a1
    Mar 27 at 7:45











  • $begingroup$
    @b3m2a1 Yeah, right. However, with ODE systems of this small dimension (dim = 2) I am not sure how to utilize vectorization since the ODE has to be solved in serial and because the vector field F may be very nonlinear.
    $endgroup$
    – Henrik Schumacher
    Mar 27 at 7:53












  • 1




    $begingroup$
    Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
    $endgroup$
    – Shinaolord
    Mar 26 at 22:07










  • $begingroup$
    You're welcome.
    $endgroup$
    – Henrik Schumacher
    Mar 26 at 22:08










  • $begingroup$
    I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
    $endgroup$
    – Shinaolord
    Mar 26 at 22:08






  • 1




    $begingroup$
    This is exactly the kind of thing I like to show people when they complain about the slowness of Mathematica. Of course with some cleverness in vectorized operations Compile could probably be avoided altogether if one so desired.
    $endgroup$
    – b3m2a1
    Mar 27 at 7:45











  • $begingroup$
    @b3m2a1 Yeah, right. However, with ODE systems of this small dimension (dim = 2) I am not sure how to utilize vectorization since the ODE has to be solved in serial and because the vector field F may be very nonlinear.
    $endgroup$
    – Henrik Schumacher
    Mar 27 at 7:53







1




1




$begingroup$
Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
$endgroup$
– Shinaolord
Mar 26 at 22:07




$begingroup$
Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
$endgroup$
– Shinaolord
Mar 26 at 22:07












$begingroup$
You're welcome.
$endgroup$
– Henrik Schumacher
Mar 26 at 22:08




$begingroup$
You're welcome.
$endgroup$
– Henrik Schumacher
Mar 26 at 22:08












$begingroup$
I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
$endgroup$
– Shinaolord
Mar 26 at 22:08




$begingroup$
I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
$endgroup$
– Shinaolord
Mar 26 at 22:08




1




1




$begingroup$
This is exactly the kind of thing I like to show people when they complain about the slowness of Mathematica. Of course with some cleverness in vectorized operations Compile could probably be avoided altogether if one so desired.
$endgroup$
– b3m2a1
Mar 27 at 7:45





$begingroup$
This is exactly the kind of thing I like to show people when they complain about the slowness of Mathematica. Of course with some cleverness in vectorized operations Compile could probably be avoided altogether if one so desired.
$endgroup$
– b3m2a1
Mar 27 at 7:45













$begingroup$
@b3m2a1 Yeah, right. However, with ODE systems of this small dimension (dim = 2) I am not sure how to utilize vectorization since the ODE has to be solved in serial and because the vector field F may be very nonlinear.
$endgroup$
– Henrik Schumacher
Mar 27 at 7:53




$begingroup$
@b3m2a1 Yeah, right. However, with ODE systems of this small dimension (dim = 2) I am not sure how to utilize vectorization since the ODE has to be solved in serial and because the vector field F may be very nonlinear.
$endgroup$
– Henrik Schumacher
Mar 27 at 7:53

















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














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





















































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







Popular posts from this blog

Adding axes to figuresAdding axes labels to LaTeX figuresLaTeX equivalent of ConTeXt buffersRotate a node but not its content: the case of the ellipse decorationHow to define the default vertical distance between nodes?TikZ scaling graphic and adjust node position and keep font sizeNumerical conditional within tikz keys?adding axes to shapesAlign axes across subfiguresAdding figures with a certain orderLine up nested tikz enviroments or how to get rid of themAdding axes labels to LaTeX figures

Tähtien Talli Jäsenet | Lähteet | NavigointivalikkoSuomen Hippos – Tähtien Talli

Do these cracks on my tires look bad? The Next CEO of Stack OverflowDry rot tire should I replace?Having to replace tiresFishtailed so easily? Bad tires? ABS?Filling the tires with something other than air, to avoid puncture hassles?Used Michelin tires safe to install?Do these tyre cracks necessitate replacement?Rumbling noise: tires or mechanicalIs it possible to fix noisy feathered tires?Are bad winter tires still better than summer tires in winter?Torque converter failure - Related to replacing only 2 tires?Why use snow tires on all 4 wheels on 2-wheel-drive cars?