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
$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 performance-tuning
$endgroup$
|
show 5 more comments
$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 performance-tuning
$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 useRule
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, 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 aTable
anyways.
$endgroup$
– Henrik Schumacher
Mar 26 at 21:40
3
$begingroup$
Why not just getNDSolve[]
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
|
show 5 more comments
$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 performance-tuning
$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
differential-equations numerical-integration performance-tuning
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 useRule
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, 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 aTable
anyways.
$endgroup$
– Henrik Schumacher
Mar 26 at 21:40
3
$begingroup$
Why not just getNDSolve[]
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
|
show 5 more comments
3
$begingroup$
AppendTo
is quadratic time complexity. Might be better to preallocate and set by index. Also it'll be much faster to not useRule
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, 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 aTable
anyways.
$endgroup$
– Henrik Schumacher
Mar 26 at 21:40
3
$begingroup$
Why not just getNDSolve[]
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
|
show 5 more comments
1 Answer
1
active
oldest
votes
$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.
$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 withCompile
, 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 operationsCompile
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 fieldF
may be very nonlinear.
$endgroup$
– Henrik Schumacher
Mar 27 at 7:53
|
show 1 more comment
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
);
);
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%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
$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.
$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 withCompile
, 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 operationsCompile
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 fieldF
may be very nonlinear.
$endgroup$
– Henrik Schumacher
Mar 27 at 7:53
|
show 1 more comment
$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.
$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 withCompile
, 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 operationsCompile
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 fieldF
may be very nonlinear.
$endgroup$
– Henrik Schumacher
Mar 27 at 7:53
|
show 1 more comment
$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.
$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.
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 withCompile
, 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 operationsCompile
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 fieldF
may be very nonlinear.
$endgroup$
– Henrik Schumacher
Mar 27 at 7:53
|
show 1 more comment
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 withCompile
, 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 operationsCompile
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 fieldF
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
|
show 1 more comment
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.
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%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');
);
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
3
$begingroup$
AppendTo
is quadratic time complexity. Might be better to preallocate and set by index. Also it'll be much faster to not useRule
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 aTable
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