0

This is a follow-up question to this one.

I am working on a project in which I have two systems of integro-differential equations. I asked this question because it seems NDSolve does not like a NIntegrate inside the equations. I tried to use the first of the solutions in this answer, but I have found a problem.

Here is part of my code for the smaller sub-problem:

G = 1;
k = 1; 
R = 1;
\[Sigma] = 1;
P[t_, r_] := R T[t, r];
e[t_, r_] := 3/2 k T[t, r] + \[Sigma] T[t, r]^4;
L[t_, r_] := -4 Pi r^2 D[T[t, r]^4, r];
eq1 = D[\[Rho][t, r], t] + D[\[Rho][t, r] v[t, r], r];
int[F_, t_?NumericQ, r_?NumericQ] := NIntegrate[F[t, r1] r1^2, {r1, 0.1, r}];
guess[t_, r_] = 1/r^2;
{dens, temp, vel} =  NestList[NDSolveValue[{D[\[Rho][t, r], t] + D[\[Rho][t, r] v[t, r], r] == 0, D[v[t, r], t] + 1/\[Rho][t, r] D[P[t, r], r] + (4 Pi G int[#, t, r])/r^2 == 0, D[e[t, r], t] + P[t, r]/\[Rho][t, r] v[t, r] (D[\[Rho][t, r], r] +  D[P[t, r], r] + \[Rho][t, r] (4 Pi G int[#, t, r])/r^2) + D[L[t, r], r] == 0, T[0, r] == 1, T[t, 0.1] == 1,  Derivative[0, 1][T][t, 0.1] == 0, v[0, r] == 1, v[t, 0.1] == 1, \[Rho][0, r] == 1, \[Rho][t, 0.1] == 1}, {\[Rho], T, v}, {t, 0, 6 10^-7}, {r, 0.1, 0.13}] &, guess,   2][[-1]]

the value of the constants is not important now, so I set them to 1. The problem is that after the first iteration, NestList tries to substitute into # the whole system of three solutions instead of just dens and obviuosly NIntegrate does not know what to do.

How do I tell NestList to substitute only the first solution?

mattiav27
  • 6,677
  • 3
  • 28
  • 64
  • @thorimur no it says: Part::partd: Part specification guess[[1]] is longer than depth of object. – mattiav27 May 19 '21 at 19:13
  • thorimur's answer is correct. A more basic fix is to modify # to #[[1]] and NestList[…, guess, 2] to NestList[…, {guess, , }, 2]. Well, reading your recent questions, I have a strong feeling that the problem you're facing isn't something can be handled with those built-in super functions only. I really suggest you to put a bit more effort to learn the core language of Mathematica. Ref: https://mathematica.stackexchange.com/q/5059/1871 – xzczd May 20 '21 at 07:31

1 Answers1

1

Using Replace[#, {x_, y__} :> x] in place of # to address the list case only seems to work:

G = 1;
k = 1;
R = 1;
\[Sigma] = 1;
P[t_, r_] := R T[t, r];
e[t_, r_] := 3/2 k T[t, r] + \[Sigma] T[t, r]^4;
L[t_, r_] := -4 Pi r^2 D[T[t, r]^4, r];
eq1 = D[\[Rho][t, r], t] + D[\[Rho][t, r] v[t, r], r];
int[F_, t_?NumericQ, r_?NumericQ] := 
  NIntegrate[F[t, r1] r1^2, {r1, 0.1, r}];
guess[t_, r_] = 1/r^2;
{dens, temp, vel} = 
 NestList[NDSolveValue[{D[\[Rho][t, r], t] + 
        D[\[Rho][t, r] v[t, r], r] == 0, 
      D[v[t, r], t] + 
        1/\[Rho][t, r] D[P[t, r], 
          r] + (4 Pi G int[Replace[#, {x_, y__} :> x], t, r])/r^2 == 
       0, D[e[t, r], t] + 
        P[t, r]/\[Rho][t, r] v[t, 
          r] (D[\[Rho][t, r], r] + 
           D[P[t, r], 
            r] + \[Rho][t, 
             r] (4 Pi G int[Replace[#, {x_, y__} :> x], t, r])/r^2) + 
        D[L[t, r], r] == 0, T[0, r] == 1, T[t, 0.1] == 1, 
      Derivative[0, 1][T][t, 0.1] == 0, v[0, r] == 1, 
      v[t, 0.1] == 1, \[Rho][0, r] == 1, \[Rho][t, 0.1] == 
       1}, {\[Rho], T, v}, {t, 0, 6 10^-7}, {r, 0.1, 0.13}] &, guess, 
   2][[-1]]

Of course if you're doing this many times, there might be more efficient ways, such as defining a function

maybefirst[x_List] := First[x]
maybefirst[x_] := x

or other similar fixes!

thorimur
  • 9,010
  • 18
  • 32