1
$\begingroup$

I have a set of four coupled Ordinary Differential Equations and I would like to solve it numerically with NDSolve. The explicit forms of these equations are not very instructive here because they are complicated (for your information, it is the field equations of a modified theory of gravity with spherical symmetry).

EDIT : Link to download my notebook : Notebook

For my purpose here, I have an initial condition at a large radius and I want to integrate backward in $r$ until $R=1.0$. Therefore, after defining my equations and my conditions, I use the following instruction

sm = NDSolve[{eqn, eqY, eqU, eqp, n[rmax] == nmax, Y[rmax] == Ymax,U[rmax] == Umax,p[rmax]==0.0}, {n, Y, U, p}, {r, 1.0, rmax},AccuracyGoal -> 20, PrecisionGoal -> 10]

The integration seems successful since I have the following output without any errors

Output of NDSolve without any errors

But then, when I want to plot the solution (for example the function $N(r)$) with the command :

plotm = Plot[Evaluate[n[r] /. sm], {r, 1.0, rmax},PlotTheme -> "Scientific", ImageSize -> Large]

there is nothing under $r\approx 5$ even though the integration domain is [1.0,rmax]. See on this picture what I get

Plot

I have the same problem by plotting $Y(r)$ or $U(r)$ (but their plots stop at different radius). It looks like the integrator failed to reach $R=1$ but usually, it prints an error message in this case. Even more weird, if I ask directly for the numerical value of $N(r)$, it seems to work perfectly. For example, with the instruction :

For[i = 0, i <= 100, i++, Block[{r}, r = 1 + 49.0/100 i; Print[r, "\t", n[r] /. sm]]]

I get the following result

Numerical values

I have tried to change the vertical axis range on the plot : it changes nothing. The plot stops below some value of r greater than 1.
I am not an expert with these features of Mathematica so may be I have forgotten something obvious but actually I do not see what it is.

Any idea ?

Thank you !

$\endgroup$
  • 2
    $\begingroup$ No full code,no idea,there will be no answer. $\endgroup$ – Mariusz Iwaniuk Mar 25 at 18:42
  • $\begingroup$ You really want all my code ? It is very long. And all the beginning is about deriving the field equations, there is no problem about this part... $\endgroup$ – Romain Gervalle Mar 25 at 18:46
  • 2
    $\begingroup$ Without code other people can run, you won't be able to get any help. $\endgroup$ – J. M.'s technical difficulties Mar 25 at 22:54
  • $\begingroup$ Ok but is it possible to put the .nb file somewhere directly ? Because I think it is more covenient than to copy-paste everything on a post. $\endgroup$ – Romain Gervalle Mar 26 at 10:00
  • $\begingroup$ I edited my post to add a link to download my notebook. $\endgroup$ – Romain Gervalle Mar 26 at 11:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.