-1
$\begingroup$

I would like to solve a set of differential equations numerically. During the process I need to use a discrete variable to annul one of my integration variables. But I also need to control what happens at r->Infinity for another variable (because for this one I have a boundary condition and not an initial condition). Mathematica says that it can not deal with discrete variables and boundary conditions. So I wanted to use a while-loop, but it doesn't work and i don't know why.

Here is my code:

β = - 4.5 ;
Γ = 2.34 ;
K = 0.0195 ;
ρ0 = 10^(18) ;
c = 3 10^8 ;
G = 6.67 10^(-11) ;
ϵ[r_] = (ρ[r] + K ρ0 (ρ[r]/ρ0)^Γ) c^2 ;
p[r_] = K ρ0 (ρ[r]/ρ0)^Γ c^2 (Γ - 1) ; 
A[r_] = Exp[(1/2) β φ[r]^2] ;    
Eqn1 = (μ'[r] - (4 Pi G/c^2) r^2 (ϵ[r]/c^2) (A[r])^4 - (1/2) r (r - 2 μ[r]) (ψ[r])^2 == 0) ;
Eqn2 = (ν'[r] - (8 Pi G/c^4) r^2 (A[r])^4 p[r]/(r - 2 μ[r]) - r (ψ[r])^2 + 2 μ[r]/(r (r - 2 μ[r])) == 0) ;
Eqn3 = (φ'[r] - ψ[r] == 0) ;    
Eqn4 = (ψ'[r] - (4 Pi G/c^4) r (A[r])^4/(r - 2 μ[r]) (β φ[r] (ϵ[r] - 3 p[r]) + r ψ[r] (ϵ[r] - p[r])) + 2 (r - μ[r]) ψ[r]/(r (r - 2 μ[r])) == 0) ;   
Eqn5 = (p'[r] + on[r]((ϵ[r] + p[r]) ((4 Pi G/c^4) r^2 (A[r])^4 p[r]/(r - 2 μ[r]) + (1/2) r (ψ[r])^2 + μ[r]/(r (r - 2 μ[r])) + β φ[r] ψ[r])) == 0) ;

a = 1 ;
n = 1 ;

While[Abs[a] > 0.0009,   
  {sol = 
    NDSolve[
      {Eqn1, Eqn2, Eqn3, Eqn4, Eqn5, 
       μ[10^(-8)] == 0, ν[10^(-8)] == 0, ρ[10^(-8)] == ρ0, ψ[10^(-8)] == 0, 
       φ[10^(-8)] == 0.5, on[10^(-8)] == 1,   
       WhenEvent[ρ[r] < 3 10^(17), on[r] -> 0],   
       WhenEvent[ρ[r] < 3*10^(17), ρ[r] -> 10^(-20)]}, 
      {μ[r], ν[r], ρ[r], ψ[r], φ[r]}, {r,10^(-8), 30000},   
      DiscreteVariables -> {on[r] ∈ {0, 1}},   
      Method -> {"StiffnessSwitching"}], a = φ[10^9] /. sol}; n++]; 
  {n, Abs[a]}  

It always gives me back {2, {0.97887}} which violates the while-loop condition.

$\endgroup$

closed as off-topic by MarcoB, Yves Klett, m_goldberg, user9660, Jason B. Mar 3 '16 at 14:42

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – MarcoB, Yves Klett, m_goldberg, Community, Jason B.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ Your Eqn4 contains a syntax error $\endgroup$ – Dr. belisarius Mar 2 '16 at 21:44
  • $\begingroup$ Could you tell me where please ? I must admit that I don't see. $\endgroup$ – arnoulxmsu Mar 2 '16 at 21:47
  • $\begingroup$ There is an unclosed parenthesis there $\endgroup$ – Dr. belisarius Mar 2 '16 at 21:48
  • $\begingroup$ You're right thanks ! but it is ok on my computer so the problem doesn't come from there $\endgroup$ – arnoulxmsu Mar 2 '16 at 21:50
  • $\begingroup$ What are you changing from one iteration to the next? n is just a counter $\endgroup$ – Dr. belisarius Mar 2 '16 at 23:05
2
$\begingroup$

This is a partial answer. I don't have time to do all your debugging. Your While-loop is syntactically unsound and has semantic problems as well.

The following edit of your loop addresses the syntax issues but is probably worthless because of the remaining semantic issues. However, it might move you forward some on your debugging.

While[Abs[a] > 0.0009,
  sol = 
    NDSolve[
     {Eqn1, Eqn2, Eqn3, Eqn4, Eqn5,
      μ[10^(-8)] == 0, ν[10^(-8)] == 0, ρ[10^(-8)] == ρ0, ψ[10^(-8)] == 0, 
      φ[10^(-8)] == 0.5, on[10^(-8)] == 1,
      WhenEvent[ρ[r] < 3 10^(17), on[r] -> 0],
      WhenEvent[ρ[r] < 3*10^(17), ρ[r] -> 10^(-20)]},
     {μ[r], ν[r], ρ[r], ψ[r], φ[r]}, {r, 10^(-8), 30000}, 
     DiscreteVariables -> {on[r] ∈ {0, 1}},
     Method -> {"StiffnessSwitching"}]; 
  a = (φ[r] /. sol /. r -> 10.^9); 
  n++]

One glaring semantic error is that you solve over the range {r, 10^(-8), 30000}, but then want to evaluate φ[r] at 10.^9, which is outside that range.

$\endgroup$

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