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
Join the discussion
Ask a follow-up here, or get your own question answered by working scientists, mathematicians and engineers — people, not an autocomplete.
Real named experts · corrections over time · the nuance an AI answer skips
11 replies · 3K views
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] == ?
 
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 a lot, 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}]
 
There are still some problems, i think there is abit of an error for some values of c(=alpha)
 
for c=1 the graph is perfect