# Finding ODE constants

I am trying to find the constant from the following ODE describing the flow in a cylindrical pipe of radius $$R$$.

$$0=-\frac{\Delta P}{L}+\mu\left(\frac{1}{r}\frac{d}{dr}\left(r\frac{d v_z}{dr} \right )\right)+\rho g \cos(\theta)$$

Where

$$\frac{d v_z}{dr}\bigg\rvert_{r=0}=0$$

And

$$v_z\bigg\rvert_{r=R}=0$$

Using the DSolve I can get the general solution.

eqn1 = {0 == -(dP/L) + \[Mu]*((1/r)*D[r*D[v[r], r], r]) + \[Rho]*g*Cos[B]};
sol1 = FullSimplify[DSolve[eqn1, v[r], r]]


We obtain the general profile with constant $$C_1$$ and $$C_2$$

{{v[r] -> C[2] + (r^2*(dP - g*L*\[Rho]*Cos[B]))/(4*L*\[Mu]) + C[1]*Log[r]}}


I can also find the specific profile by including the boundary conditions

eqn1 = {0 == -(dP/L) + \[Mu]*((1/r)*D[r*D[v[r], r], r]) + \[Rho]*g*Cos[\[Beta]], v[R] == 0, (D[v[r], r] /. r -> 0) == 0};
sol1 = FullSimplify[DSolve[eqn1, v[r], r]]


To obtain the resolved velocity profile

{{v[r] -> ((r - R) (r + R) (dP -
g L \[Rho] Cos[\[Beta]]))/(4 L \[Mu])}}


How can I obtain the values for constant $$C_1$$ and $$C_2$$ using Mathematica.

• You have already done so when imposing the boundary conditions. It's not clear what you are asking
– bmf
Commented Jan 22, 2023 at 5:09
• btw, your latex notation does not agree with your code. I used your code. In the code you have v[R] == 0 but in Latex you wrote $v_z\bigg\rvert_{r=R}=0$ which looks like v'[R]==0 and not v[R]==0 Commented Jan 22, 2023 at 5:26

eqn1 = 0 == -(z/L) + μ*((1/r)*D[r*D[v[r], r], r]) + ρ*g*Cos[B]
ic = {v[R] == 0, v'[0] == 0}
sol1 = DSolveValue[eqn1, v[r], r]


Take derivative

D[sol1, r]


For bounded solution at $$r=0$$, C[1] has to be zero (since one of your initial conditions is v'[0] == 0. Hence the solution now becomes

sol2 = sol1 /. C[1] -> 0


To find C[2] apply the second initial conditions

eq = (sol2 /. r -> R) == 0
c2 = Solve[eq, C[2]][[1, 1]]


These are your C[1] and C[2].

To verify, we can plugin in these in the general solution, and compare with the solution given by Mathematica

eqn1 = 0 == -(z/L) + μ*((1/r)*D[r*D[v[r], r], r]) + ρ*g*Cos[B]
ic = {v[R] == 0, v'[0] == 0}
mmaSolution = DSolveValue[{eqn1, ic}, v[r], r]


sol1 = sol1 /. {C[1] -> 0, c2}


mmaSolution == sol1 // Simplify


• Than you so much Nasser. This is exactly what I needed. Commented Jan 23, 2023 at 1:09