Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Help in Mathematica code for solutions expansion of differential equations

  1. Oct 16, 2011 #1
    About serval differential equations where A, B, D, g, \chi, c are functions of r
    &-\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&
    The solution near r=0 is given by
    & 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\\
    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.
  2. jcsd
  3. Oct 16, 2011 #2
    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:

    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.
    Last edited: Oct 16, 2011
  4. Oct 16, 2011 #3
    Thanks a lot.
    I will try and report it later.
  5. Oct 17, 2011 #4
    May I ask a stupid question?
    Are d0 and d1 functions of rstart or just a number using rstart=0.001?

  6. Oct 18, 2011 #5
    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}]


    Attached Files:

  7. Oct 19, 2011 #6
    First problem in your equations.
    r^2 (c[r])^2 g[r]) a[r] = 0
    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
    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.
    Last edited: Oct 19, 2011
  8. Oct 19, 2011 #7
    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:

    Code (Text):

    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}]
  9. Oct 20, 2011 #8
    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.
    Last edited: Oct 20, 2011
  10. Oct 26, 2011 #9
    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?

    Attached Files:

    Last edited: Oct 26, 2011
  11. Oct 27, 2011 #10
    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.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook