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

Time travel short story where a man arrives in the late 19th century in a time machine and then sends the machine back into the past

Opposite of a diet

Can I Retrieve Email Addresses from BCC?

Do I need a multiple entry visa for a trip UK -> Sweden -> UK?

Applicability of Single Responsibility Principle

What is difference between behavior and behaviour

How do I define a right arrow with bar in LaTeX?

Greatest common substring

Should my PhD thesis be submitted under my legal name?

Tiptoe or tiphoof? Adjusting words to better fit fantasy races

Best way to store options for panels

Can I use my Chinese passport to enter China after I acquired another citizenship?

Finding all intervals that match predicate in vector

Ways to speed up user implemented RK4

How to be diplomatic in refusing to write code that breaches the privacy of our users

Are there any comparative studies done between Ashtavakra Gita and Buddhim?

Is the destination of a commercial flight important for the pilot?

voltage of sounds of mp3files

Is there a problem with hiding "forgot password" until it's needed?

How can a jailer prevent the Forge Cleric's Artisan's Blessing from being used?

Modify casing of marked letters

How does residential electricity work?

How do I keep an essay about "feeling flat" from feeling flat?

What to do with wrong results in talks?



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!










share|improve this question









$endgroup$







  • 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
    2 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
    2 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
    2 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
    2 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
    2 hours ago















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!










share|improve this question









$endgroup$







  • 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
    2 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
    2 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
    2 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
    2 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
    2 hours ago













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!










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






share|improve this question













share|improve this question











share|improve this question




share|improve this question










asked 3 hours ago









ShinaolordShinaolord

808




808







  • 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
    2 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
    2 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
    2 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
    2 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
    2 hours ago












  • 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
    2 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
    2 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
    2 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
    2 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
    2 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
2 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
2 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
2 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
2 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
2 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
2 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
2 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
2 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
2 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
2 hours ago










1 Answer
1






active

oldest

votes


















4












$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







share|improve this answer









$endgroup$












  • $begingroup$
    Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
    $endgroup$
    – Shinaolord
    2 hours ago










  • $begingroup$
    You're welcome.
    $endgroup$
    – Henrik Schumacher
    2 hours ago










  • $begingroup$
    I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
    $endgroup$
    – Shinaolord
    2 hours ago










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









4












$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







share|improve this answer









$endgroup$












  • $begingroup$
    Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
    $endgroup$
    – Shinaolord
    2 hours ago










  • $begingroup$
    You're welcome.
    $endgroup$
    – Henrik Schumacher
    2 hours ago










  • $begingroup$
    I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
    $endgroup$
    – Shinaolord
    2 hours ago















4












$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







share|improve this answer









$endgroup$












  • $begingroup$
    Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
    $endgroup$
    – Shinaolord
    2 hours ago










  • $begingroup$
    You're welcome.
    $endgroup$
    – Henrik Schumacher
    2 hours ago










  • $begingroup$
    I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
    $endgroup$
    – Shinaolord
    2 hours ago













4












4








4





$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







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








share|improve this answer












share|improve this answer



share|improve this answer










answered 2 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
    2 hours ago










  • $begingroup$
    You're welcome.
    $endgroup$
    – Henrik Schumacher
    2 hours ago










  • $begingroup$
    I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
    $endgroup$
    – Shinaolord
    2 hours ago
















  • $begingroup$
    Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
    $endgroup$
    – Shinaolord
    2 hours ago










  • $begingroup$
    You're welcome.
    $endgroup$
    – Henrik Schumacher
    2 hours ago










  • $begingroup$
    I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
    $endgroup$
    – Shinaolord
    2 hours ago















$begingroup$
Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
$endgroup$
– Shinaolord
2 hours ago




$begingroup$
Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!
$endgroup$
– Shinaolord
2 hours ago












$begingroup$
You're welcome.
$endgroup$
– Henrik Schumacher
2 hours ago




$begingroup$
You're welcome.
$endgroup$
– Henrik Schumacher
2 hours ago












$begingroup$
I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
$endgroup$
– Shinaolord
2 hours ago




$begingroup$
I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.
$endgroup$
– Shinaolord
2 hours ago

















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

Are there any AGPL-style licences that require source code modifications to be public? Planned maintenance scheduled April 23, 2019 at 23:30 UTC (7:30pm US/Eastern) Announcing the arrival of Valued Associate #679: Cesar Manara Unicorn Meta Zoo #1: Why another podcast?Force derivative works to be publicAre there any GPL like licenses for Apple App Store?Do you violate the GPL if you provide source code that cannot be compiled?GPL - is it distribution to use libraries in an appliance loaned to customers?Distributing App for free which uses GPL'ed codeModifications of server software under GPL, with web/CLI interfaceDoes using an AGPLv3-licensed library prevent me from dual-licensing my own source code?Can I publish only select code under GPLv3 from a private project?Is there published precedent regarding the scope of covered work that uses AGPL software?If MIT licensed code links to GPL licensed code what should be the license of the resulting binary program?If I use a public API endpoint that has its source code licensed under AGPL in my app, do I need to disclose my source?

2013 GY136 Descoberta | Órbita | Referências Menu de navegação«List Of Centaurs and Scattered-Disk Objects»«List of Known Trans-Neptunian Objects»

Mortes em março de 2019 Referências Menu de navegação«Zhores Alferov, Nobel de Física bielorrusso, morre aos 88 anos - Ciência»«Fallece Rafael Torija, o bispo emérito de Ciudad Real»«Peter Hurford dies at 88»«Keith Flint, vocalista do The Prodigy, morre aos 49 anos»«Luke Perry, ator de 'Barrados no baile' e 'Riverdale', morre aos 52 anos»«Former Rangers and Scotland captain Eric Caldow dies, aged 84»«Morreu, aos 61 anos, a antiga lenda do wrestling King Kong Bundy»«Fallece el actor y director teatral Abraham Stavans»«In Memoriam Guillaume Faye»«Sidney Sheinberg, a Force Behind Universal and Spielberg, Is Dead at 84»«Carmine Persico, Colombo Crime Family Boss, Is Dead at 85»«Dirigent Michael Gielen gestorben»«Ciclista tricampeã mundial e prata na Rio 2016 é encontrada morta em casa aos 23 anos»«Pagan Community Notes: Raven Grimassi dies, Indianapolis pop-up event cancelled, Circle Sanctuary announces new podcast, and more!»«Hal Blaine, Wrecking Crew Drummer, Dies at 90»«Morre Coutinho, que editou dupla lendária com Pelé no Santos»«Cantor Demétrius, ídolo da Jovem Guarda, morre em SP»«Ex-presidente do Vasco, Eurico Miranda morre no Rio de Janeiro»«Bronze no Mundial de basquete de 1971, Laís Elena morre aos 76 anos»«Diretor de Corridas da F1, Charlie Whiting morre aos 66 anos às vésperas do GP da Austrália»«Morreu o cardeal Danneels, da Bélgica»«Morreu o cartoonista Augusto Cid»«Morreu a atriz Maria Isabel de Lizandra, de "Vale Tudo" e novelas da Tupi»«WS Merwin, prize-winning poet of nature, dies at 91»«Atriz Márcia Real morre em São Paulo aos 88 anos»«Mauritanie: décès de l'ancien président Mohamed Mahmoud ould Louly»«Morreu Dick Dale, o rei da surf guitar e de "Pulp Fiction"»«Falleció Víctor Genes»«João Carlos Marinho, autor de 'O Gênio do Crime', morre em SP»«Legendary Horror Director and SFX Artist John Carl Buechler Dies at 66»«Morre em Salvador a religiosa Makota Valdina»«مرگ بازیکن‌ سابق نساجی بر اثر سقوط سنگ در مازندران»«Domingos Oliveira morre no Rio»«Morre Airton Ravagniani, ex-São Paulo, Fla, Vasco, Grêmio e Sport - Notícias»«Morre o escritor Flavio Moreira da Costa»«Larry Cohen, Writer-Director of 'It's Alive' and 'Hell Up in Harlem,' Dies at 77»«Scott Walker, experimental singer-songwriter, dead at 76»«Joseph Pilato, Day of the Dead Star and Horror Favorite, Dies at 70»«Sheffield United set to pay tribute to legendary goalkeeper Ted Burgin who has died at 91»«Morre Rafael Henzel, sobrevivente de acidente aéreo da Chapecoense»«Morre Valery Bykovsky, um dos primeiros cosmonautas da União Soviética»«Agnès Varda, cineasta da Nouvelle Vague, morre aos 90 anos»«Agnès Varda, cineasta francesa, morre aos 90 anos»«Tania Mallet, James Bond Actress and Helen Mirren's Cousin, Dies at 77»e