Help with troubleshooting code for system of PDEs

Question

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?

Answer 1

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

Answer 2

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.