# NDSolve is slow when defined as a function

I am solving a set of coupled PDEs, dependent variables (n,l) which are initial value problems on the domain x: {-10,10} and time: {0,1}. Let's say we have an external function u and some initial data:

u[x_] := Sin[2 Pi/20 x];
n00[x_] := 1;
l00[x_] := 0;


I can type out the equations directly:

{ns1, ls1} = NDSolve[{
D[n[x, t], t] - D[n[x, t], x, x] +
D[1/2 u[x] n[x, t] + 1/4 D[u[x], x] l[x, t], x] - 1 + n[x, t] ==
NeumannValue[0, x == -10] + NeumannValue[0, x == 10],

D[l[x, t], t] - D[l[x, t], x, x] + u[x] n[x, t] + 2 l[x, t] +
1/2 u[x] D[l[x, t], x] ==
NeumannValue[0, x == -10] + NeumannValue[0, x == 10],

n[x, 0] == n00[x],
l[x, 0] == l00[x]
}
, {n, l}, {t, 0, 1}, {x, -10, 10}][[1]] // AbsoluteTiming


This gives me a time of approximately 0.05 seconds. However, if I calculate the same equations but define the NDSolve as a function:

Propagatenl[tstart_, tend_, ninit_, linit_] := NDSolve[{
D[n[x, t], t] - D[n[x, t], x, x] +
D[1/2 u[x] n[x, t] + 1/4 D[u[x], x] l[x, t], x] - 1 + n[x, t] ==
NeumannValue[0, x == -10] + NeumannValue[0, x == 10],

D[l[x, t], t] - D[l[x, t], x, x] + u[x] n[x, t] + 2 l[x, t] +
1/2 u[x] D[l[x, t], x] ==
NeumannValue[0, x == -10] + NeumannValue[0, x == 10],

n[x, tstart] == ninit[x],
l[x, tend] == linit[x]
}
, {n, l}, {t, tstart, tend}, {x, -10, 10}][[1]];


And then say

{ns1, ls1} = Propagatenl[0, 1, n00, l00] // AbsoluteTiming


This will give me a time of 0.13 seconds, more than twice as long. Is there a reason why? Is there a way to use NDSolve in a function without increasing the calculation time? Thank you.

• it looks like just the overhead of building the definition of the function Propagatenl first time it is used. Because if you call this function again and again, now it has same timing as the other one because the function is now defined. Btw, your use of the letter l for function names and as first letter makes your code very hard to read as it looks like 1. Jan 9, 2023 at 23:47
• I don't think so? I've done: {ns1, ls1} = Propagatenl[0, 1, n00, l00] // AbsoluteTiming {ns1, ls1} = Propagatenl[0, 1, n00, l00] // AbsoluteTiming. And both times, it gives around 0.1 seconds.
– xdu
Jan 10, 2023 at 16:57

I don't think so? I've done: {ns1, ls1} = Propagatenl[0, 1, n00, l00] // AbsoluteTiming {ns1, ls1} = Propagatenl[0, 1, n00, l00] // AbsoluteTiming. And both times, it gives around 0.1 seconds.

Too large for comment. Here is what I get on V 13.2. I put screen shots as it is easier to show (click on the screen shot to enlarge). I am using same code you posted.

Let compare the above without using the function

## Summary

When using a function, first time there is more overhead (0.23 vs. 0.19), but after that, the overhead is much less, since it is defined and lookup is fast now. It takes around 0.07-0.08 seconds now to make the calls vs. 0.23 for first time. With no function it takes about 0.05 seconds. So there is still very small difference due to the function overhead.

Make sure to do this each time from a clean kernel as shown above to compare.

This is all on V 13.2 on windows 10