3
$\begingroup$

i'm currently trying to numerically solve a set of coupled ODE's to obtain the functions p(r), h(r) and m(r) in the range of r1 <= r <= r2 with initial conditions m(r1)=a=const and p(r1)=b=const. It is my goal to vary the parameters a and b to find the highest possible value for m(r2). So far i have no problems solving the ODE's for fixed a and b, but i have no idea how to vary the two parameters. I'm really new to Mathematica and after working around with ParametricNDSolve, ParametricPlot etc. for one week without making any progress i am now desperate for help.

Here's the code i made so far that solves the ODE's where i fixed r1=1.8, r2=2.2 and a=0.8, b=1. Please note that there is also another condition that a < r1/2, but i don't think that this is very important to my problem.

r1 = 1.8;
r2 = 2.2;
ode1 = {h'[r] == (1/r)*(1 - 8*Pi*p[r]*r^2 - h[r])};
ode2 = {p'[r] == -(2 p[r]/r) - (h'[r]/h[r])*p[r]*((1 + 8*Pi*p[r]*r^2 - 3*h[r])/(1 - 8*Pi*p[r]*r^2 - h[r]))};
{H, P, M} = {h, p, m} /.NDSolve[{ode1, ode2 , m[r] == (r/2)*(1 - h[r]), m[r1] == 0.8, p[r1] == 1}, {h, p, m}, {r, r1, r2}] // First

Thank you very much in advance.

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – user9660 May 3 '16 at 12:05
2
$\begingroup$

Oh! I'm coming to late!

r1 = 1.8;
r2 = 2.2;
ode1 = h'[r] == (1/r)*(1 - 8*Pi*p[r]*r^2 - h[r]);
ode2 = p'[
    r] == -(2 p[r]/r) - (h'[r]/h[r])*
     p[r]*((1 + 8*Pi*p[r]*r^2 - 3*h[r])/(1 - 8*Pi*p[r]*r^2 - h[r]));

sol = ParametricNDSolve[{ode1, ode2, m[r] == (r/2)*(1 - h[r]), 
   m[r1] == a, p[r1] == b}, {h, p, m}, {r, r1, r2}, {a, b}]

enter image description here

One way to find max m[a,b][r2]

MaximalBy[
 Table[Evaluate[{a, b, m[a, b][2.2] /. sol}], {a, 0.2, 0.8, 0.1}, {b, 
    0, 0.5, 0.1}][[1]], Last]
(* {{0.2, 0.3, 0.970011}} *)

Manipulate[
 Plot[Evaluate[{h[a, b][r], p[a, b][r], m[a, b][r]} /. sol], {r, r1, 
   r2}, PlotTheme -> "Detailed", PlotLegends -> {h, p, m}], {a, 0.1, 
  0.3, Appearance -> "Open"}, {b, 0.1, 0.5, Appearance -> "Open"}]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you very, very much for your answer! Don't worry, even though you answered later, your answer sure provides more useful information to me since you also gave me the means to find the maximum value! I really appreciate your help, you sure did me an enormous favor! Even though i need to do some research to understand what exactly you did, but reverse-engineering won't be too hard i guess :p $\endgroup$ – Xeno May 3 '16 at 17:17
  • $\begingroup$ @Xeno Of course, you can do it $\endgroup$ – user36273 May 3 '16 at 20:12
3
$\begingroup$
r1 = 1.8;
r2 = 2.2;
ode1 = {h'[r] == (1/r)*(1 - 8*Pi*p[r]*r^2 - h[r])};
ode2 = {p'[r] == -(2 p[r]/r) - (h'[r]/h[r])*
      p[r]*((1 + 8*Pi*p[r]*r^2 - 3*h[r])/(1 - 8*Pi*p[r]*r^2 - h[r]))};

hpm = ParametricNDSolveValue[{ode1, ode2, m[r] == (r/2)*(1 - h[r]), 
   m[r1] == a, p[r1] == b}, {h, p, m}, {r, r1, r2}, {a, b}]

Plot[Evaluate[Through[hpm[0.8, 1.][t]]], {t, 1.8, 2.2}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Thank you very, very much for your help! Through, ParametricNDSolveValue? Those certainly are new commands to me! Thank you again for your help, i really appreciate it. $\endgroup$ – Xeno May 3 '16 at 17:09

Your Answer

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

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