0
$\begingroup$
NDSolve[
  {D[r[dm], dm] == 1/(4π r[dm]^2 ρ2[dm]), 
   D[P[dm], dm] == (-G dm)/(4π r[dm]^4),
   D[L[dm], dm] == ϵ,
   P[dm] == (ρ2[dm]/(μ mu ) kb T2[dm])+(a/3 T2[dm]^4),
   D[T2[dm], dm] == (3 κ L[dm])/(256 π^2 σ r[dm]^4 T2[dm]^3)},
  {r[dm], L[dm], T2[dm], ρ2[dm], P[dm]}, {dm, Mstart, 0.01}]

I want to solve this system of coupled equations, but as you can see it is without boundary conditions. Therefore, when I try to add one like:

L[Mstart] = L0

like so

NDSolve[
  {D[r[dm], dm] == 1/(4π r[dm]^2 ρ2[dm]), 
   D[P[dm], dm] == (-G dm)/(4π r[dm]^4),
   D[L[dm], dm] == ϵ,
   P[dm] == (ρ2[dm]/(μ mu ) kb T2[dm])+(a/3 T2[dm]^4),
   D[T2[dm], dm] == (3 κ L[dm])/(256 π^2 σ r[dm]^4 T2[dm]^3)},
   L[Mstart] = L0,
  {r[dm], L[dm], T2[dm], ρ2[dm], P[dm]}, {dm, Mstart, 0.01}]

This is still not very well defined, but I get the following error message: (Note that Mstart and L0 are real numbers)

NDSolve::dsfun: L$4740[1.5912*10^30] == 3.846*10^26 cannot be used as a function.

The problem with this is that in examples on the reference site they do do it this way, but with me it doesn't seem to work, I'm doing this within a module, so perhaps that has something to do with it?

Edit

The parameters are all numbers,

G = 6.674*10^-11;
kb = 1.38064852*10^-23;
σ = 5.670367*10^-8;
mu = 1.66053904*10^-27;
μ = 1.28617;
ϵ = 400;
a=7.56572*10^-16

κ can be set to 1.

Edit

The starting values, or boundary conditions are as follows:

L[Mstart]==L0,
r[Mstart]==Rstart,
\[Rho]2[Mstart]==\[Rho]start,
T2[Mstart]==Tstart

With

L0=1.5912*10^30
Rstart=5.0112*10^8 
\[Rho]=7180.8
Tstart=5.7*10^6

Update

When I use the method kindly provided by zhk, but with the desired boundary conditions(As given in the second Edit) I recieve the following error:

NDSolve::ndsz: At dm$7416 == 1.5879824736744075`*^30, step size is effectively zero; singularity or stiff system suspected.

When I added a dynamic [Kappa] it seemed to work a little, but after realizing a mistake, it just gives the same error, so I doubt it has something to do with [Kappa], but for reference, here is the [Kappa] I used:

data1=Import["http://www.uio.no/studier/emner/matnat/astro/AST3310/v18/beskjeder/opacity.txt","Table"];
data2=Drop[data1,{2}];
datalogR=Drop[data2[[1]],1];
datalogT=Drop[data2[[All,1]],1];
datalog\[Kappa]=Drop[Drop[data2,1],0,1];
\[Kappa]function=ListInterpolation[datalog\[Kappa],{{3.75,8},{-8,1}}];

\[Kappa][dm_]=(10^\[Kappa]function[Log[10,T2[dm]],Log[10,\[Rho]2[dm]/(T2[dm]/10^6)^3]])/10

The correction was that I forgot the cube added to the (T2[dm]/10^6) part, perhaps the increased value on [Kappa] was making a possible step size greater then it already is?

$\endgroup$
  • $\begingroup$ I would be good if you can add the missing parameters, otherwise people will most likely not want to try. $\endgroup$ – user21 Mar 1 '18 at 14:05
  • 1
    $\begingroup$ The problem I see is that you gave the boundary condition as separate argument when it should be included in the 1st argument, the list of equations. $\endgroup$ – m_goldberg Mar 1 '18 at 15:02
1
$\begingroup$
G = 6.674*10^-11;
kb = 1.38064852*10^-23;
\[Sigma] = 5.670367*10^-8;
mu = 1.66053904*10^-27;
\[Mu] = 1.28617;
\[Epsilon] = 400;
Mstart = 0;
L0 = 1;
a = 1;
\[Kappa] = 1;

P[dm_] = (\[Rho]2[dm]/(\[Mu] mu) kb T2[dm]) + (a/3 T2[dm]^4)

Eq1 = D[r[dm], dm] == 1/(4 \[Pi] r[dm]^2 \[Rho]2[dm])

Eq2 = D[P[dm], dm] == (-G dm)/(4 \[Pi] r[dm]^4)

Eq3 = D[L[dm], dm] == \[Epsilon]

Eq4 = D[T2[dm], dm] == (3 \[Kappa] L[dm])/(256 \[Pi]^2 \[Sigma] r[dm]^4 T2[dm]^3)

sol = NDSolve[{Eq1, Eq2, Eq3, Eq4, r[Mstart] == 1, L[Mstart] == L0, 
   T2[Mstart] == 1, \[Rho]2[Mstart] == 1}, {r[dm], L[dm], 
   T2[dm], \[Rho]2[dm]}, {dm, Mstart, 0.01}]

Plot[{r[dm], L[dm], T2[dm], \[Rho]2[dm]} /. sol, {dm, Mstart, 0.01}]

enter image description here

$\endgroup$
  • $\begingroup$ Well, at least I know now that it can work, which is a relief, and I forgot to give values for all the "start" values, and the boundary conditions, so I will update my post with them. $\endgroup$ – Hidde Rinsema Mar 2 '18 at 9:41

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.