# Help with troubleshooting code for system of PDEs

I'm trying to solve the following system of equations

$$\frac{\partial \rho}{\partial t} + \frac{f(t)}{r^{1/2}}\frac{\partial \rho}{\partial r} - \frac{1}{2} \rho \frac{f(t)}{r^{3/2}} = 0$$

$$b\rho = \left[2\frac{df}{dt}r^{1/2}\right]^2 - \left(\frac{f(t)}{r^{1/2}} \right)^2$$

I want to use NDSolve to solve for $$\rho(r,t)$$ and $$f(t)$$ numerically, subject to the conditions: $$f(1)=1$$, $$\rho(r,0) = \frac{1}{\sqrt{2\pi}\sigma}\text{exp}\left(-\frac{(r-R)^2}{2\sigma^2}\right)$$ and that $$\rho(r,t)$$ goes to zero as $$r\rightarrow \infty$$ (I put in $$\rho(100 R,t)$$. The mathematica code looks like the one below:

R = 100;
b = 0.81;
\[Sigma] = 1;

eq1 = D[\[Rho][r, t], t] + f[t]/r^(1/2) D[\[Rho][r, t], r] -
1/2 \[Rho][r, t] f[t]/r^(3/2) == 0;
eq2 = b \[Rho][r, t] == (2 D[f[t], t] r^(1/2))^2 - (f[t]/r^(1/2))^2;

bcs = {\[Rho][r, 0] ==
1/((2 \[Pi])^(1/2) \[Sigma])
Exp[-((r - R)^2/(2 \[Sigma]^2))]; \[Rho][10 R, t] == 0;
f[1] == 1};

vars = {\[Rho][r, t], f[t]};

sol = NDSolve[{eq1, eq2, bcs}, vars, {r, 10, 10 R}, {t, 0, 10 R}]



When I run it, Mathematica doesnt solve it and instead throws out the following errors/warnings:

Function::fpct: Too many parameters in {r,t} to be filled from Function[{r,t},1][t].

NDSolve::overdet: There are fewer dependent variables, {[Rho][r,t]}, than equations, so the system is overdetermined.

How do I fix this?

• For NDSolve to solve a system of PDE, all functions have to be declared as dependent on all variables, that may be expressed here by $$\{ ..., \partial_r f [r,t]==0\}, }{\rho[r,t ], f[r, t]\}$$ Commented Oct 5, 2023 at 18:36
• In bcs you need to use a comma instead of a semicolons. Commented Oct 5, 2023 at 19:07
• @RolandF where do I put \partial_r f [r,t]==0? Under the boundary conditions or as a third equation?
– jboy
Commented Oct 5, 2023 at 19:30

Extended comment:

R = 100;
b = 0.81;
\[Sigma] = 1;

eq1 = D[\[Rho][r, t], t] + f[t]/r^(1/2) D[\[Rho][r, t], r] -1/2 \[Rho][r, t] f[t]/r^(3/2) == 0;
eq2 = b \[Rho][r, t] == (2 D[f[t], t] r^(1/2))^2 - (f[t]/r^(1/2))^2;


From eq2 we get \[Rho][r, t]

rho = Function[{r, t}, Values@Solve[eq2, \[Rho] [r, t]][[1, 1]] // Evaluate]
(*Function[{r, t}, 1.23457 (-(f[t]^2/r) + 4. r Derivative[1][f][t]^2)]*)


Substituting this result we get an ode for f[t]

odef = eq3[[1]] /. \[Rho] -> rho // Simplify
(*(0.1875 f[t]^3)/Sqrt[r] +
r f[t] Derivative[1][f][t] (-0.25 + 0.25 Sqrt[r] Derivative[1][f][t]) +
1. r^3 Derivative[1][f][t] (f^\[Prime]\[Prime])[t] == 0*)


This ode depends on r,t which contradicts the assumption "f[t] depends only on t"!

• I understand what you're saying. My first attempt was similar; I rewrote the replacement rule to replace $\rho$ and it's derivatives, and was left with an equation in terms of f[t] and r. But it wasn't clear how to proceed. Choose value of r individually and generate a family of solutions? It seemed to require that the problem be reformulated, and maybe that's the correct way.
– user87932
Commented Oct 6, 2023 at 16:26
• The problem is, what happens to the boundary conditions, two of which are defined in terms of $\rho$? The replacement plugs an expression with f and f' in, so it's not clear how to untangle it and generate a usable b.c.
– user87932
Commented Oct 6, 2023 at 16:40
• @jdp You could try to solve eq1 for f[t] and substitude this result into eq2. This gives a pde in \[Rho][r, t] . Perhaps now it's easier to handle the boundary conditions? Commented Oct 6, 2023 at 16:57
• That's a good idea. I won't have time till later today to try it. Though the 3rd equation in the bcs is f[1]==1, so that could still be a potential problem.
– user87932
Commented Oct 6, 2023 at 18:00

EDIT: Since this answer kept growing in length after multiple edits, I went back and deleted stuff that's no longer relevant, though I did retain the explanations of what I found and did.

As Daniel Huber pointed out, you have a typo in bcs. I noticed that your eq2 specifies $$\rho(r,t)$$ as a function of f[t] and r, so I used Solve to get an expression for $$\rho(r,t)$$ and applied the solution to eq1, replacing the $$\rho(r,t)$$ term with the solution from eq2 in terms of f,f',r - but leaving the $$\rho$$ derivatives untouched, as well as the boundary conditions.

That took care of the first problem but now I ran into a new one - the exponential in your boundary condition has a huge value for large r, and you get a numerical underflow.

While thinking about this, I looked at a suggestion made by Ulrich Neumann to try solving eq1 for f[t] in terms of rho and replacing f[t],f'[t] everywhere. Sort of the inverse of what I did previously. The problem is that the boundary conditions get hopelessly tangled and I see no way to translate the result into the relatively simple f(1) = 1 specified. So I resumed trying to make progress with what I had originally.

The comment under the question about f,rho requiring the same arguments is correct; the easy way to see this is to try DSolve without boundary conditions. It complains and points to the appropriate documentation page to explain the message. I rewrote the equation and added a term to have the r-derivative of f be 0. I also split the "R" into r0 (where to center the Gaussian initially), and rtmax - the maximum value for the t,r grid. I also renamed the variables sigma and rho because I got weary of hitting escape to type Greek letters, and they don't come out well when I try to paste code here anyway, and made some formatting changes to reduce the number of parenthesis.

I made some changes to your initial condition for rho. One earlier problem was that you had inconsistent boundary conditions, which caused NDSolve to bail and return unevaluated. You specified rho[r,0] to be a Gaussian, but also stipulate that rho[rtMax,t] == 0. What happens at {r,t} = {rtMax,0}? Both conditions apply at that point, and they're not consistent; the Gaussian is small, but not zero. If you have an inconsistency in your bc's at even one point, NDSolve issues a message and exits. I wrapped your Gaussian in an If statement and force it to zero beyond some distance from r0. This may lead to a small discontinuity, but it shouldn't have too much effect.

I did get past a few of the earlier problems only to run into a new set, generally suggesting further problems with boundary conditions. After further thought, it occurred to me that it's unrealistic to say that f[t] (now f[r,t]) remains set indefinitely, and that perhaps a useful boundary condition to add is to arrange for it to turn off somewhere at the edge of the grid. So I added this as another condition, along with If statements to modify the f,f derivative boundary conditions to be consistent; i.e. start dropping to the value at the region edge and to avoid conflicts where conditions overlap in r,t:

With[{r0 = 10, b = .81, sigma = 1, rtMax = 100},
eq = -.61728395061722839 f[r,
t] (-f[r, t]^2/r + 4 r Derivative[0, 1][f][r, t]^2)/Sqrt[r^3] +
Derivative[0, 1][rho][r, t] +
f[r, t] Derivative[1, 0][rho][r, t]/Sqrt[r] == 0;
bcs = {rho[r, 0] ==
If[Abs[r - r0] < 10 sigma,
1/(Sqrt[2 Pi] sigma) Exp[-((r - r0)^2/(2 sigma^2))], 0],
rho[rtMax, t] == 0, f[r, 0] == 1, f[rtMax, t] == If[t > 0, .9, 1],
Derivative[1, 0][f][r, t] == If[r < (rtMax - 1), 0, -.1]};
sol = NDSolve[{eq, bcs}, {rho, f}, {r, r0, rtMax}, {t, 0, rtMax}]]


I now get a solution for both rho and f over the full grid with no error messages. Here's the output and a couple plots:

I was trying a number of things to get this/anything to work, so you may find that some of the changes I made aren't compatible with your original problem, or wish to make different choices for how to handle the behavior of f along the r-axis, etc. But at least you have a working baseline now to use as a reference point for further modifications.

• Hi, so I was thinking maybe it can be handled by adjusting the values of the parameter sigma to be larger... I havent tried it yet because I'm having trouble replicating the results when I try out your code. Specifically, in bcs /. r-> 10R, when I do this it just throws out {f[1.]==1}, the first two conditions regarding \rho are not returned
– jboy
Commented Oct 6, 2023 at 3:51
• @jdp Interesting approach, but the substitution deriving eq3 shouldn't contain \[Rho][r, t] anymore! I get (0.1875 f[t]^3)/Sqrt[r] + r f[t] Derivative[1][f][ t] (-0.25 + 0.25 Sqrt[r] Derivative[1][f][t]) + 1. r^3 Derivative[1][f][t] (f^\[Prime]\[Prime])[t] == 0 for eq3, which unfortunately depends on t and r. Commented Oct 6, 2023 at 8:02
• @jdp That means (eq3 depends on f,r,t), assumption f[t] seems to be wrong! Commented Oct 6, 2023 at 10:27
• @UlrichNeumann I just reran everything with a fresh kernel and I got the same expression for eq3 which I posted above, which looks quite different from yours. I'm not sure what's going on. What I did was to eliminate the $\rho(r,t)$ term from eq1 and replaced it with f[t],f'[t], and r terms, but left all the $\rho$ derivatives alone. Whether this is the correct way to handle this problem is unclear, but it seems to get past the OP's first problem at the bottom of the question, before running into an underflow in the boundary conditions.
– user87932
Commented Oct 6, 2023 at 15:55
• @jboy I just reran with a fresh kernel and reproduced what I posted above; I don't know why you're having trouble replicating it. If you only get the third term, the first two disappeared - but off to where? What if you just type bcs? Do you see all three terms there? I tried bumping sigma up to 10^5, but still have problems. Something keeps generating non-numerical derivative errors at the upper bound.
– user87932
Commented Oct 6, 2023 at 16:05