Ways to speed up user implemented RK4

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!

asked 4 hours ago

Shinaolord

858

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
4 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
4 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
4 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
4 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
4 hours ago

|
show 1 more comment

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!

asked 4 hours ago

Shinaolord

858

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
4 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
4 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
4 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
4 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
4 hours ago

|
show 1 more comment

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!

asked 4 hours ago

Shinaolord

858

rk4[f_, variables_, valtinit_, tinit_, tfinal_, nsteps_] := 

  Module[{table, xlist, ylist, step, k1, k2, k3, k4}, 

   xlist = tinit;

   step = N[(tfinal - tinit)/(nsteps)];

   ylist = valtinit;

   table = {{xlist, ylist}};

   Table[

    k1 = step* f /. MapThread[Rule, {variables, ylist}]; (* 

    Equivalent to step* f/.Thread[Rule[variables,ylist]]*)

    k2 = step*f /. MapThread[Rule, {variables, k1/2 + ylist}];

    k3 = step*f /. MapThread[Rule, {variables, k2/2 + ylist}];

    k4 = step*f /. MapThread[Rule, {variables, k3 + ylist}];

    ylist += 1/6 (k1 + 2 (k2 + k3) + k4);

    xlist += step;

    AppendTo[table, {xlist, ylist}];

    {xlist, ylist}, nsteps];

   table

   ];

Example Input:

funclist = {-x + y, x - y};

initials = {1, 2};

variables = {x, y};

init = 0;

final = 200;

nstep = 20000;

approx = rk4[funclist, variables, initials, init, final, nstep]//AbsoluteTiming;

{3.59932,{...}}

I'd love some suggestions!

differential-equations numerical-integration

asked 4 hours ago

Shinaolord

858

asked 4 hours ago

Shinaolord

858

asked 4 hours ago

Shinaolord

858

asked 4 hours ago

Shinaolord

858

asked 4 hours ago

Shinaolord

858

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
4 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
4 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
4 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
4 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
4 hours ago

|
show 1 more comment

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
4 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
4 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
4 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
4 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
4 hours ago

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.

– b3m2a1
4 hours ago

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!

– Shinaolord
4 hours ago

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.

– Henrik Schumacher
4 hours ago

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.

– Henrik Schumacher
4 hours ago

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

– Shinaolord
4 hours ago

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

– Shinaolord
4 hours ago

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

– Shinaolord
4 hours ago

|
show 1 more comment

1 Answer
1

active

oldest

votes

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

For given stepsize τ and given vectorfield F, this creates a CompiledFunction cStep that computes a single Runge-Kutta step

F = X [Function] {-Indexed[X, 2], Indexed[X, 1]};



τ = 0.01;

Block[{YY, Y, k1, k2, k3, k4},



  YY = Table[Compile`GetElement[Y, i], {i, 1, 2}];

  k1 = τ F[YY];

  k2 = τ F[0.5 k1 + YY];

  k3 = τ F[0.5 k2 + YY];

  k4 = τ F[k3 + YY];



  cStep = With[{code = YY + (k1 + 2. (k2 + k3) + k4)/6. },

    Compile[{{Y, _Real, 1}},

     code,

     CompilationTarget -> "C",

     RuntimeOptions -> "Speed"

     ]

    ]

  ];

Now we can apply it 2 million times with NestList and still need only 2 seconds.

nsteps = 20000000;

xlist = Range[0., step nsteps, step];

Ylist = NestList[cStep, initials, nsteps]; // AbsoluteTiming // First

2.08678

answered 4 hours ago

Henrik Schumacher

57.9k579159

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

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

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

add a comment |

Your Answer

StackExchange.ifUsing("editor", function () {
return StackExchange.using("mathjaxEditing", function () {
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\$","\$"]]);
});
});
}, "mathjax-editing");

StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "387"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});

function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: false,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});

}
});

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

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

Name

Required, but never shown

1 Answer
1

active

oldest

votes

1 Answer
1

active

oldest

votes

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

For given stepsize τ and given vectorfield F, this creates a CompiledFunction cStep that computes a single Runge-Kutta step

F = X [Function] {-Indexed[X, 2], Indexed[X, 1]};



τ = 0.01;

Block[{YY, Y, k1, k2, k3, k4},



  YY = Table[Compile`GetElement[Y, i], {i, 1, 2}];

  k1 = τ F[YY];

  k2 = τ F[0.5 k1 + YY];

  k3 = τ F[0.5 k2 + YY];

  k4 = τ F[k3 + YY];



  cStep = With[{code = YY + (k1 + 2. (k2 + k3) + k4)/6. },

    Compile[{{Y, _Real, 1}},

     code,

     CompilationTarget -> "C",

     RuntimeOptions -> "Speed"

     ]

    ]

  ];

Now we can apply it 2 million times with NestList and still need only 2 seconds.

nsteps = 20000000;

xlist = Range[0., step nsteps, step];

Ylist = NestList[cStep, initials, nsteps]; // AbsoluteTiming // First

2.08678

answered 4 hours ago

Henrik Schumacher

57.9k579159

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

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

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

add a comment |

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

For given stepsize τ and given vectorfield F, this creates a CompiledFunction cStep that computes a single Runge-Kutta step

F = X [Function] {-Indexed[X, 2], Indexed[X, 1]};



τ = 0.01;

Block[{YY, Y, k1, k2, k3, k4},



  YY = Table[Compile`GetElement[Y, i], {i, 1, 2}];

  k1 = τ F[YY];

  k2 = τ F[0.5 k1 + YY];

  k3 = τ F[0.5 k2 + YY];

  k4 = τ F[k3 + YY];



  cStep = With[{code = YY + (k1 + 2. (k2 + k3) + k4)/6. },

    Compile[{{Y, _Real, 1}},

     code,

     CompilationTarget -> "C",

     RuntimeOptions -> "Speed"

     ]

    ]

  ];

Now we can apply it 2 million times with NestList and still need only 2 seconds.

nsteps = 20000000;

xlist = Range[0., step nsteps, step];

Ylist = NestList[cStep, initials, nsteps]; // AbsoluteTiming // First

2.08678

answered 4 hours ago

Henrik Schumacher

57.9k579159

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

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

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

add a comment |

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

For given stepsize τ and given vectorfield F, this creates a CompiledFunction cStep that computes a single Runge-Kutta step

F = X [Function] {-Indexed[X, 2], Indexed[X, 1]};



τ = 0.01;

Block[{YY, Y, k1, k2, k3, k4},



  YY = Table[Compile`GetElement[Y, i], {i, 1, 2}];

  k1 = τ F[YY];

  k2 = τ F[0.5 k1 + YY];

  k3 = τ F[0.5 k2 + YY];

  k4 = τ F[k3 + YY];



  cStep = With[{code = YY + (k1 + 2. (k2 + k3) + k4)/6. },

    Compile[{{Y, _Real, 1}},

     code,

     CompilationTarget -> "C",

     RuntimeOptions -> "Speed"

     ]

    ]

  ];

Now we can apply it 2 million times with NestList and still need only 2 seconds.

nsteps = 20000000;

xlist = Range[0., step nsteps, step];

Ylist = NestList[cStep, initials, nsteps]; // AbsoluteTiming // First

2.08678

answered 4 hours ago

Henrik Schumacher

57.9k579159

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

For given stepsize τ and given vectorfield F, this creates a CompiledFunction cStep that computes a single Runge-Kutta step

F = X [Function] {-Indexed[X, 2], Indexed[X, 1]};



τ = 0.01;

Block[{YY, Y, k1, k2, k3, k4},



  YY = Table[Compile`GetElement[Y, i], {i, 1, 2}];

  k1 = τ F[YY];

  k2 = τ F[0.5 k1 + YY];

  k3 = τ F[0.5 k2 + YY];

  k4 = τ F[k3 + YY];



  cStep = With[{code = YY + (k1 + 2. (k2 + k3) + k4)/6. },

    Compile[{{Y, _Real, 1}},

     code,

     CompilationTarget -> "C",

     RuntimeOptions -> "Speed"

     ]

    ]

  ];

Now we can apply it 2 million times with NestList and still need only 2 seconds.

nsteps = 20000000;

xlist = Range[0., step nsteps, step];

Ylist = NestList[cStep, initials, nsteps]; // AbsoluteTiming // First

2.08678

answered 4 hours ago

Henrik Schumacher

57.9k579159

answered 4 hours ago

Henrik Schumacher

57.9k579159

answered 4 hours ago

Henrik Schumacher

57.9k579159

answered 4 hours ago

Henrik Schumacher

57.9k579159

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

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

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

add a comment |

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

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

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

Damn, you definitely know how to use Mathematica A LOT more efficiently than I do. Thanks!

– Shinaolord
4 hours ago

You're welcome.

– Henrik Schumacher
4 hours ago

I'll have to play around with Compile, it definitely seems like a massive speed up if used correctly.

– Shinaolord
4 hours ago

add a comment |

draft saved

draft discarded

Thanks for contributing an answer to Mathematica Stack Exchange!

Please be sure to answer the question. Provide details and share your research!

But avoid …

Asking for help, clarification, or responding to other answers.

Making statements based on opinion; back them up with references or personal experience.

Use MathJax to format equations. MathJax reference.

To learn more, see our tips on writing great answers.

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

Post as a guest

Name

Required, but never shown

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

Name

Required, but never shown

Name

Required, but never shown

This page is only for reference, If you need detailed information, please check here

搜尋此網誌

Gfrktyl