Help in Mathematica code for solutions expansion of differential equations

[jadyliber]

About serval differential equations where A, B, D, g, \chi, c are functions of r
\begin{eqnarray}
&-\frac{{\chi}'}{r}+\frac{c'}{c}\left(\frac{g'}{g} -{\chi}'\right)=\frac{e^{\chi}(q A B)^2}{r^2 g^2 c^2}& \\
&c c''+c c'\left(\frac{g'}{g}+\frac{2}{r} -\frac{{\chi}'}{2} \right)=-\frac{B'^2}{2 r^2}+\frac{e^{\chi}}{2 g^2 r^2}(q A B)^2 & \\
&-g'\left(\frac{1}{r}+\frac{c'}{2c}\right)- g \left(\frac{1}{r^2}+\frac{3c'}{c r}+\frac{c''}{c} \right)+3 = \frac{e^{\chi}}{4}A'^2+\frac{g B'^2}{4 r^2 c^2}+\frac{e^{\chi} (q A B)^2}{4 g r^2 c^2}+\frac{e^{\chi}}{4 }D'^2 & \\
&A'' +\left(\frac{2}{r}+\frac{\chi'}{2}+\frac{c'}{c} \right) A'-\frac{q^2B^2}{r^2 c^2 g } A=0& \\
&B'' +\left(\frac{g'}{g}-\frac{\chi'}{2}-\frac{c'}{c} \right) B'+\frac{e^{\chi}q^2A^2}{g } B=0& \\
&D'' +\left(\frac{2}{r}+\frac{\chi'}{2}+\frac{c'}{c} \right) D'=0&
\end{eqnarray}
The solution near r=0 is given by
\begin{eqnarray}
& A\sim A_0 e^{-\alpha/r}, ~~~B\sim B_0\left(1-\frac{e^{\chi_0}q^2A_0^2}{4\alpha^2}e^{-2\alpha/r}\right),
~~~c\sim c_0\left(1+\frac{e^{\chi_0}A_0^2}{8r^2}e^{-2\alpha/r}\right), & \nonumber\\
& \chi\sim \chi_0-\frac{e^{\chi_0}A_0^2\alpha}{2r^3}e^{-2\alpha/r},~~~g\sim r^2-\frac{e^{\chi_0}A_0^2\alpha}{2r}e^{-2\alpha/r}, ~~~D\sim -\frac{e^{\chi_0}A_0^2\alpha}{8r^3}e^{-2\alpha/r}. & \nonumber\\
\end{eqnarray}
where terms with subscript 0 are constant.

My question is how to get the whole range numerical solution of r using mathematica. Any code for similar equations are welcome.

 

[jackmell]

You just need to plug them into NDSolve, assign all the constants and initial conditions, and then just turn the crank. However, don't use upper-case letters for user-defined variables. For example "D" is a build-in function in Mathematica. Just use d[r] or myd[r]. I mean you can use all those greek letters but that just takes more typing. I'll do the one for d[r]. Note the use of double equal signs and I'm starting close to r=0:

rstart=0.001;
mysol=NDSolve[{d''[r]+(2/r+x'[r]/2+c'[r]/c[r]) d'[r]==0, xequation, cequation, gequation, aequation, bequation,d[rstart]==d0,d'[rstart]==d1, all the other initial conditions here},{d,x,c,g,a,b},{r,rstart,rend}]

However, don't expect to get it the first time. May have to tweak the parameters to get it to run without encounting precision problems, singular problems, or other numerical limitations. I suggest for starters, just let all the constants equal to one.

 

[jadyliber]

Thanks a lot.
I will try and report it later.

 

[jadyliber]

May I ask a stupid question?
Are d0 and d1 functions of rstart or just a number using rstart=0.001?
Thanks!

You just need to plug them into NDSolve, assign all the constants and initial conditions, and then just turn the crank. However, don't use upper-case letters for user-defined variables. For example "D" is a build-in function in Mathematica. Just use d[r] or myd[r]. I mean you can use all those greek letters but that just takes more typing. I'll do the one for d[r]. Note the use of double equal signs and I'm starting close to r=0:

rstart=0.001;
mysol=NDSolve[{d''[r]+(2/r+x'[r]/2+c'[r]/c[r]) d'[r]==0, xequation, cequation, gequation, aequation, bequation,d[rstart]==d0,d'[rstart]==d1, all the other initial conditions here},{d,x,c,g,a,b},{r,rstart,rend}]

However, don't expect to get it the first time. May have to tweak the parameters to get it to run without encounting precision problems, singular problems, or other numerical limitations. I suggest for starters, just let all the constants equal to one.

 

[jadyliber]

The following is the code, and it does not work. I do not know whether the code itself is ok. Thanks
"
rstart = 0.001;
mysol = NDSolve[{d''[r] + (2/r + x'[r]/2 + c'[r]/c[r]) d'[r] ==
0, -(x'[r]/r) + c'[r]/c[r] (g'[r]/g[r] - x'[r]) == (
q^2 (b[r])^2 (a[r])^2)/(r^2 (c[r])^2 g[r]),
c[r] c''[r] +
c[r] c'[r] (g'[r]/g[r] + 2/r - x'[r]/2) == -((b'[r])^2/(
2 r^2 )) + ( e^x[r] q^2 (b[r])^2 (a[r])^2)/(
2 r^2 (g[r])^2), -g'[r] (1/r + c'[r]/(2 c[r] )) -
g[r] (2/r^2 + (3 c'[r])/(r c[r]) - c''[r]/c[r]) +
3 == -(( e^x[r] (a'[r])^2)/4) + (g[r] (b'[r])^2)/(
4 r^2 (c[r])^2) ( e^x[r] q^2 (b[r])^2 (a[r])^2)/(
4 r^2 g[r] (c[r])^2) + ( e^x[r] (d'[r])^2)/4,
a''[r] + (2/r + x'[r]/2 + c'[r]/c[r]) a'[r] - ( q^2 (b[r])^2)/(
r^2 (c[r])^2 g[r]) a[r] = 0,
b''[r] + (g'[r]/g[r] - x'[r]/2 - c'[r]/c[r]) b'[r] + (
e^x[r] q^2 (a[r])^2)/g[r] b[r] = 0,
d[rstart] == -(q/(8 (rstart)^3)) e^(-((2 q)/rstart)),
d'[rstart] == -(q^2/(4 (rstart)^5)) e^(-((2 q)/rstart)),
x[rstart] == -(q/(2 (rstart)^3)) e^(-((2 q)/rstart)),
x'[rstart] == -(q^2/(rstart)^5) e^(-((2 q)/rstart)),
c[rstart] == 1, c'[rstart] == q/(4 (rstart)^4) e^(-((2 q)/rstart)),
g[rstart] == rstart^2, g'[rstart] == 2 rstart,
a[rstart] == e^(-(q/rstart)),
d'[rstart] == q/(rstart)^2 e^(-(q/rstart)), b[rstart] == 1,
b'[rstart] == q/(2 (rstart)^2) e^(-((2 q)/rstart)) }, {d, x, c, g,
a, b}, {r, rstart, rend}]
"

You just need to plug them into NDSolve, assign all the constants and initial conditions, and then just turn the crank. However, don't use upper-case letters for user-defined variables. For example "D" is a build-in function in Mathematica. Just use d[r] or myd[r]. I mean you can use all those greek letters but that just takes more typing. I'll do the one for d[r]. Note the use of double equal signs and I'm starting close to r=0:

rstart=0.001;
mysol=NDSolve[{d''[r]+(2/r+x'[r]/2+c'[r]/c[r]) d'[r]==0, xequation, cequation, gequation, aequation, bequation,d[rstart]==d0,d'[rstart]==d1, all the other initial conditions here},{d,x,c,g,a,b},{r,rstart,rend}]

However, don't expect to get it the first time. May have to tweak the parameters to get it to run without encounting precision problems, singular problems, or other numerical limitations. I suggest for starters, just let all the constants equal to one.

 

[Bill Simpson]

First problem in your equations.
r^2 (c[r])^2 g[r]) a[r] = 0
and
e^x[r] q^2 (a[r])^2)/g[r] b[r] = 0
both need == instead of =.

Second problem.
I'm guessing, since you don't include e in your list of variables to solve for, that you mean Euler's constant when you write e.
Mathematica requres you write that as E instead of e.

Third problem.
Even when I make those fixes it complains that you have more initial conditions (12) than the order of your system (10) and it bails out.

Fourth problem.
When I inspect your notebook I find two different d`[rstart]==...
Did you possibly mean one of those to be a`[rstart]==...?
Even if I try making that change it still complains about more initial conditions than the order of your system.

Fifth problem.
NDSolve is likely not going to be happy with your q that doesn't appear to be assigned any constant value. Similarly for your rend. Typically NDSolve requires everything but the things you are solving for to be constants or have been assigned constant numeric values.

At this point I think I've made enough guesses about changes that I will stop and let you go through this, correct any errors and then try to decide how to deal with more initial conditions than the order of your system.

As a completely meaningless experiment, I inserted
rend=0.1;q=1.0;
commented out x`[rstart]==...
commented out g`[rstart]==...
and almost instantly got back 6 interpolation functions as a solution to your system.
My choice of values and selection of initial conditions to discard are almost certainly incorrect.
But that at least gives some hope you might get a solution if you use sensible changes.

 

[jackmell]

Thanks for that Bill. Sorry I left jady with a bear. Been busy. Looks like Jady needs to maybe practice with just a single DE or maybe two to get the hang of entering the correct syntax into Mathematica. You need to enter the DEs and all the initial conditions. If you have a second-order DE, then you need two initial conditions, y(ystart) and y'(ystart). If it's just a first order, then only need y(start) for the initial condtions. And as Bill mentioned, for numerical results, you have to specify values for all the constants. Also, the e thing too. You could just specify Exp[x(r)] and I noticed some typos from the equations you posted and the code you specified. I cleaned up the code, just let q=-0.1 for starters, and assigned random initial conditioons for all the initial conditions. I gave c(r) a relatively large value since it's in the denominator. Now, just supply the actual initial condtiions you wish to run the code:

rstart = 0.001; 
rend = 1; 
q = -0.1; 

mysol = NDSolve[{Derivative[2][d][r] + (2/r + Derivative[1][x][r]/2 + Derivative[1][c][r]/c[r])*
       Derivative[1][d][r] == 0, -(Derivative[1][x][r]/r) + (Derivative[1][c][r]/c[r])*
       (Derivative[1][g][r]/g[r] - Derivative[1][x][r]) == E^x[r]*((q*b[r]*a[r])^2/(r*c[r]*g[r])^2), 
    c[r]*Derivative[2][c][r] + c[r]*Derivative[1][c][r]*(Derivative[1][g][r]/g[r] + 2/r - Derivative[1][x][r]/2) == 
     -Derivative[1][b][r]^2/(2*r^2) + (E^x[r]*q*b[r]*a[r])^2/(2*r*g[r])^2, 

    (-Derivative[1][g][r])*(1/r + Derivative[1][c][r]/(2*c[r])) - 
      g[r]*(1/r^2 + (3*Derivative[1][c][r])/(r*c[r]) + Derivative[2][c][r]/c[r]) + 3 == 
     (Exp[x[r]]/4)*Derivative[1][a][r]^2 + (g[r]*Derivative[1][b][r]^2)/(4*r^2*c[r]^2) + 
      (Exp[x[r]]*(q*a[r]*b[r])^2)/(4*g[r]*r^2*c[r]^2) + (Exp[x[r]]/4)*Derivative[1][d][r]^2, 

    Derivative[2][a][r] + (2/r + Derivative[1][x][r]/2 + Derivative[1][c][r]/c[r])*Derivative[1][a][r] - 
      ((q^2*b[r]^2)/(r^2*c[r]^2*g[r]))*a[r] == 0, 

    Derivative[2][b][r] + (Derivative[1][g][r]/g[r] - Derivative[1][x][r]/2 - Derivative[1][c][r]/c[r])*
       Derivative[1][b][r] + ((Exp[x[r]]*q^2*a[r]^2)/g[r])*b[r] == 0, 

d[rstart] == 0, 
    Derivative[1][d][rstart] == 0.1,

 x[rstart] == 0.1,

 c[rstart] == 3.1,
 Derivative[1][c][rstart] == 0.2, 

    g[rstart] == 1.5,

 a[rstart] == 0,
 Derivative[1][a][rstart] == 0.1,

 b[rstart] == 2.1, 
    Derivative[1][b][rstart] == 0.1},

 {d, x, c, g, a, b}, {r, rstart, rend}];

Plot[{d[r], x[r], c[r], g[r], a[r], b[r]} /. mysol, {r, rstart, rend}]

 

[jadyliber]

Thanks for Bill's and Jack's comments and sorry for my mistake in nb first. Next I will try and report it later. So waiting for me. Thanks again.
By the way, to Bill, all your guesses are absolutely right.

 

[jadyliber]

Thanks for Jack'S and Bill's kind help. I can run some basic code mow.
While I still have some related problems.

Since I have solution near the r-> 0^+, and I still have some symmetry to set A_0=1, chi_0=0 and c_0=1, the only changing initial condition B_0 and q. That is in the following code, for given rstart and q, the only unfixed parameter is b0. In my problem, I have some other condition that b=0 as r goes to infinity. Thus for given q and b0 , in order to set b=0 at infinity r, there should be a relation between q and b0. How to find it ?

For each times given q and rstart, I should calculate d[rstart] , d'[rstart],
x[rstart], c'[rstart], a[rstart], a'[rstart],
b'[rstart] independently using initial solution and some derivatives. Is there any method to pass it? (Attachment2)

BTW, how to express nb in the reply rather than attachment?

 

[jackmell]

To post Mathematica code, first select it, then choose Cell/Convert To/Raw Input form. That changes everything to ASCII which you can then cut and paste into the forum and then choose the "Wrap code tags" from the menu bar.