I am trying to solve a complicated system of integro-differential equations, but a MWE of my problem is the following:

fun[rho_?NumericQ, r_?NumericQ] := r NIntegrate[rho r1^2, {r1, 0, r}];
NDSolve[{R'[r] == D[fun[R[r], r], r], R[0] == 1}, R, {r, 0, 1}]

The code runs forever without producing any output: how can I solve this kind of equations?

EDIT a solution for this equation is to derive once more in $r$ and substitute back R'[r], but my system is much more complicated and I do not know the initial values for the derivatives.

EDIT2 I have tried the following code:

NDSolve[{R'[r]==D[r NIntegrate[r1^2 R[r1],{r1,0,r}],r],R[0]==1},R,{r,0,1}]

but Mathematica complaints that r us not a valid number for an extreme of integration.

How can I solve theis problem?

  • $\begingroup$ It took less than 3 minutes on my Windows machine in v12.2 to produce the resulting interpolating function. $\endgroup$
    – flinty
    Commented May 16, 2021 at 15:21
  • $\begingroup$ For the toy example we can even funa[rho_, r_] = r Integrate[rho r1^2, {r1, 0, r}]; DSolve[{R'[r] == D[funa[R[r], r], r], R[0] == 1}, R, r], but this isn't suitable for the real case, right? $\endgroup$
    – xzczd
    Commented May 16, 2021 at 15:22
  • $\begingroup$ @xzczd No the real case is more complicated $\endgroup$
    – mattiav27
    Commented May 16, 2021 at 15:23
  • $\begingroup$ @flinty Just tested in v9.0.1, it takes 38 seconds to compute, and the result is correct. (Oh I forgot Derivative can handle things like Derivative[1, 0][fun][1., 1.] to some degree, related: mathematica.stackexchange.com/q/196998/1871 ) mattiav27, what version are you in? $\endgroup$
    – xzczd
    Commented May 16, 2021 at 15:27
  • $\begingroup$ @xzczd Now that I think about it: if I use fun[R[r], r], doesNIntegrate understand that R[r] must be evaluated at the point r1 inside the integral? $\endgroup$
    – mattiav27
    Commented May 16, 2021 at 15:30

3 Answers 3


If the integral is assumed to be inrooted, I can think out 2 approaches. First one is the relaxation method:

intfunc[R_, r_?NumericQ] := 
  NIntegrate[r1^2*R[r1], {r1, 0, r}, Method -> {Automatic, SymbolicProcessing -> 0}];

sollst = NestList[
    NDSolveValue[{R'[r] == intfunc[#, r] + r^3 #[r], R[0] == 1}, R, {r, 0, 2}] &, 
    Identity, 10]; // AbsoluteTiming
(* {0.402924, Null} *)

aplot = Plot[(E^(r^4/4) (-r Gamma[-(1/4)] + (r^4)^(1/4) Gamma[-(1/4), r^4/4]))/(
  4 Sqrt[2]), {r, 0, 2}, PlotRange -> All];

ListPlot[sollst[[-1]], PlotRange -> All, PlotStyle -> Red]~Show~aplot

Mathematica graphics

The next one is finite difference method (FDM), which is harder to code but more efficient at least in this case. I'll use pdetoae to facilatate the generation of difference equations:

int[expr_, {var_, 0, rlst_List}] := int[expr, {var, 0, #}] & /@ rlst
int[expr_, {var_, 0, r_?NumericQ}] := 
 With[{pos = Rescale[r, domain, {1, points}]}, 
  trap[Function @@ {var, expr}, grid[[;; pos]]]]

(* trapezoidal rule: *)
trap[f_, grid_] := 
  With[{h = -Subtract @@ domain/(points - 1)}, 
   h (Total[f /@ grid] - 1/2 (f@grid[[1]] + f@grid[[-1]]))];

eq = R'[r] == D[r intmid[r], r];
ic = R[0] == 1;

points = 80;
domain = {0, 2};
grid = Array[# &, points, domain];
difforder = 4;

(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
ptoafunc = pdetoae[{R, intmid}[r], grid, difforder];

ae = Rest@ptoafunc[eq] /. Thread[ptoafunc[intmid[r]] -> int[r1^2 R[r1], {r1, 0, grid}]];

guess[r_] = 1;
nsol = ListInterpolation[
    FindRoot[Flatten@{ae, ic}, 
      Table[{R@r, guess[r]}, {r, grid}]][[All, -1]], {grid}]; // AbsoluteTiming
(* {0.103995, Null} *)

ListPlot[nsol, PlotRange -> All, PlotStyle -> Red]~Show~aplot

Mathematica graphics

  • $\begingroup$ How should I modify the first code if R was a function of two variables? $\endgroup$
    – mattiav27
    Commented May 17, 2021 at 15:14
  • $\begingroup$ @mattiav27 That depends. Can you show a specific example? $\endgroup$
    – xzczd
    Commented May 18, 2021 at 3:11
  • $\begingroup$ An example closer to my case could be D[R[t,r,x],t]+D[R[t,r,x],r]+D[R[t,r,x],x]==D[r NIntegrate[r1^2 R[t,r1,x1],{r1,0,r},{x1,0,1}] ,r] with IC: R[t,0,x]==1 and R[0,r,x]==1. t is time. The real case is a system of 5 coupled PDE. $\endgroup$
    – mattiav27
    Commented May 18, 2021 at 5:41
  • $\begingroup$ I forgot an initial condition R[t,r,0]==1 $\endgroup$
    – mattiav27
    Commented May 18, 2021 at 5:47
  • $\begingroup$ @mattiav27 Then just define 2 helper functions. (Observe the output of D[r Integrate[r1^2 R[t, r1, x1], {r1, 0, r}, {x1, 0, 1}], r] and think about what should be defined. ) Also, notice the calculation will be much slower for PDE system. $\endgroup$
    – xzczd
    Commented May 18, 2021 at 7:10

TL;DR I'm not sure if either of my attempts are correct, will update later.

Attempt 1:

First SE post I've written so I apologize in advance if anything is out of order. The error messages from the OP coming from NIntegrate could be because NIntegrate usually has a hard time dealing with symbolic limits of integration, so NIntegrate was being called before a value was plugged into $r$.

I was able to get a solution to the problem in MM version 12.2 with:

    fun[rho_, r_] := Times[r, Integrate[Times[rho, p, p], {p, 0, r}]];
    sol = NDSolve[{R'[r] == D[fun[R[r], r], r], R[0] == 1}, R, {r, 0, 1}][[1,1,2]];
    Plot[sol[r], {r, 0, 1}]

gives a plot that looks like:Solution to an integral equation

I would note that the solution seems to be faster without using NumericQ on the function, but depending on your application, you may or may not wish to keep it there.

I was able to get two different analytical solutions, the first being:

    a1 = DSolve[{R'[r] == D[r * Integrate[q * q * R[q], {q, 0, r}], r], R[0] == p}, R[r], r]

and a more refined version over the unit interval:

    a2 = DSolve[{R'[r] == D[r * Integrate[q * q * R[q], {q, 0, r}], r], R[0] == p}, R[r], {r, 0, 1}]

yielding the expressions: \begin{equation} R\left( r\right)=\frac{1}{4\sqrt{2}} e^{\frac{r^{4}}{4}}p\left[ r\;\Gamma\left( -\frac{1}{4}\right) -\left( r^{4}\right)^{\frac{1}{4}} \Gamma\left( -\frac{1}{4}, \frac{r^{4}}{4}\right) \right] \end{equation}

where $\Gamma\left(\cdot\right)$ and $\Gamma\left(\cdot,\cdot\right)$ are intended to mean the gamma and upper-incomplete gamma functions respectively, and on the unit interval:

\begin{equation} R\left(r\right)=\frac{1}{4}e^{\frac{r^{4}}{4}}p\left[E_{5/4}\left(\frac{r^{4}}{4}\right) + 2\sqrt{2}\;r\;\Gamma\left(\frac{3}{4}\right)\right] \end{equation}

where $E_{n}$ is the exponential integral function.

To be true(er) to the original problem though, I was able to use:

    mfun[rho_?NumericQ, r_?NumericQ] := r NIntegrate[rho r1^2, {r1,0,r}];
    msol = NDSolve[{R'[r] == D[mfun[R[r], r], r], R[0] == 1}, R, {r,0,1}][[1,1,2]]

which when plotted gives the same visual as shown above.

  • $\begingroup$ I'm sorry, but as mentioned in the comment above, the NDSolve solution is incorrect, you're calculating Integrate[q * q * R[r], {q, 0, r}] rather than Integrate[q * q * R[q], {q, 0, r}] in this case, and if you Show the plots together, the difference between sol and a1/a2 is obvious. $\endgroup$
    – xzczd
    Commented May 17, 2021 at 5:04
  • $\begingroup$ I derived a second solution, will try to add in the AM. $\endgroup$ Commented May 17, 2021 at 7:08
  • $\begingroup$ @NathanWhite Very interesting answer! I'm struggeling to derive a pure numerical answer using NDSolve&NIntegrate but didn't succeed. Any idea? Thanks! $\endgroup$ Commented May 17, 2021 at 9:52


If I understand the question right, OP is looking for a numerical solution where integral isn't known analytically.

First we need a function definition of the integral part, where a function ip is passed as argument.

int[ip_][r_?NumericQ] :=Module[{u}, NIntegrate[u^2 ip[u], {u, 0, r}]  ]

Differentiating the integralequation "manually" (don't know why NDSolve can't handle D[r int[Function[{r}, R[r]]][r ],r] ) gives the ode

RR = NDSolveValue[{R'[r] == int[Function[{r}, R[r]]][r ] + r^3 R[r],R[0] == 1}, R, {r, 0, 2}]
Plot[RR[r], {r, 0, 1}]

enter image description here

The solution agrees quite well (also in the extended domain 0<r<2 , see @xzczd comment, Thanks) with the analytical one in @Nathan White`s answer:

sol=DSolve[{R'[r] == D[r * Integrate[q * q * R[q], {q, 0, r}], r], R[0] == 1}, R[r],  r ]
Plot[{RR[r], R[r] /. sol}, {r, 0, 2}, AxesLabel -> {"r","R[r]"},PlotRange -> All, PlotStyle -> {Blue, {Dashed, Red}}]

enter image description here

  • $\begingroup$ I'm sorry, but this answer is wrong. NDSolve has actually parsed your ODE to something amount to R'[r] == Integrate[u^2 R[r], {u, 0, r}] + r^3 R[r]. And if you enlarge the domain to {r, 0, 2} in NDSolve, you'll see it's still different from the analytic solution. $\endgroup$
    – xzczd
    Commented May 17, 2021 at 12:50
  • $\begingroup$ @xzczd Astonishing, thanks. I only checked the range {r,0, 1}. That means, NDSolve parses the wrong ode. I have already noticed that the function call only works in the form int[Function[{r}, R[r]]][r ] .Variants int[Function[{u}, R[u]]][r ] or int[R][r] or int[R[#]&][r] fail . $\endgroup$ Commented May 17, 2021 at 13:10
  • $\begingroup$ @xzczd I modified the definition int[ip_][r_?NumericQ] :=Module[{u}, NIntegrate[u^2 ip[u], {u, 0, r}] ]. Solution now also agrees with the analytical one in the range {r,0,2}! $\endgroup$ Commented May 17, 2021 at 13:20
  • $\begingroup$ Well, did you forgot to modify the domain of NDSolve or Plot? The solution is still wrong according to my test… $\endgroup$
    – xzczd
    Commented May 17, 2021 at 13:43
  • $\begingroup$ No I modified the domain. I 'll add the comparison. $\endgroup$ Commented May 17, 2021 at 13:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.