How can I plot b[a]/b[0.0001] vs. a using Mathematica?

  • Thread starter Thread starter kptsilva
  • Start date Start date
  • Tags Tags
    Mathematica Plot
kptsilva
Messages
33
Reaction score
0
hey guys,
I came across a research paper stating to numerically integrate the following equation.

2/3 a^2 b''[a] + (1 - w[a]) a b'[a] - (1 + w[a]) (1 - 3 c w[a]) b[a]

A Boundary condition is given b'[0.0001]=0
Where;

w[a_] := 2*a^(3*(1 + c))/(1 + 2*a^(3*(1 + c)));

(c=1) (c is a variable but let's consider a particular instance c=1)

'a' goes from 10^-4 to 1000 in a log scale.

I want to plot b[a]/b[0.0001] vs. a

how can i plot this using mathematica?
 
Physics news on Phys.org
2/3 a^2 b''[a] + (1 - w[a]) a b'[a] - (1 + w[a]) (1 - 3 c w[a]) b[a] == ?
 
?==0
 
See the attached notebook. The second plot on a log scale is probably more instructive.
 

Attachments

the paper says that they have numerically integrated that equation.I can't understand that statement. The program should be written with NIntegrate rather than NDSolve right?
 
Since it's a second order equation, you would need to integrate it twice. Integrating along step-by-step is really what NDSolve is doing. I don't see a simple way to write a solution using the NIntegarate command. Does it matter?
 
Figures 2 and 3 are the same graph, just plotted on different scales. They've numerically solved the differential equation just like I did. Also, in your original post, you forgot the minus sign in w[a]. This makes a big difference! Now it looks like the paper. See attached.
 

Attachments

phyzguy thanks alot, in my program i have forgotten the negative sign it gives me the graphs perfectly.
thhis is the code i wrote,
c = 1;
w[a_] := (2*a^(3*(1 + c)))/(1 + (2*a^(3*(1 + c))));
fun = (2/3)*(a^2)*b''[a] + (1 - w[a])*a*
b'[a] - (1 + w[a])*(1 - 3*c*w[a])*b[a]; sol =
DSolve [{fun == 0, b'[10^-4] == 0}, b, a];
A = LogLinearPlot[Evaluate[b[a]/b[10^-4] /. sol], {a, 10^-4, 10^3},
PlotRange -> {-500, 12000}]
 
  • #10
There are still some problems, i think there is abit of an error for some values of c(=alpha)
 
  • #11
for c=1 the graph is perfect
 
  • #12
Any news to my problem?
 
Back
Top