I am trying to solve this differential equation for a heat transfer problem: \begin{equation} kt\frac{\partial^2 T}{\partial x^2} = \epsilon \sigma T^4, \ \ \ T(0) = T_0, \ \ \ \frac{\partial T}{\partial x} \Big|_{x=L} = 0 \end{equation}
where $k$, $t$, $\epsilon$, $\sigma$, $T_0$ and $L$ are constants.
Mathematica's NDSolve won't touch this. However, substituting with the dimensionless $\Theta = \frac{T}{T_0}$ and $\rho = \frac{x}{L}$, the problem becomes
\begin{equation} A\frac{\partial^2 \Theta}{\partial \rho^2} = \Theta^4, \ \ \ \Theta(0) = 1, \ \ \ \frac{\partial \Theta}{\partial \rho} \Big|_{\rho=1} = 0 \end{equation}
Where $A = \frac{kt}{L^2 T_0^3 \epsilon \sigma}$ is a constant (feel free to double check this, but I am pretty confident with the math).
Pluging this into Mathematica like so
T0 = 238;
L = 1.5;
A = 0.001356;
tau = NDSolve[{A*Th''[r] == Th[r]^4, Th'[1] == 0, Th[0] == 1}, Th, {r, 0, 1}];
T[x_] = Evaluate[T0*Th[x/L]/.tau];
Plot[T[x] - 273, {x, 0, L}, PlotRange -> All, AxesLabel -> {"x [m]", "T [°C]"}]
gives me
NDSolve::ndsz: At r == 0.0530353097865862`, step size is effectively zero; singularity or stiff system suspected. >>
However, with the wrong equation (setting $A=1$):
tau = NDSolve[{Th''[r] == Th[r]^4, Th'[1] == 0, Th[0] == 1}, Th, {r, 0, 1}]
I get a nice, wrong temperature profile.
What gives? I feel like this is a math issue rather than a coding one, but I may be wrong.
I there a mathematical route to using the $\frac{\partial^2 \Theta}{\partial \rho^2} = \Theta^4$ solution and scale it with $A$ somehow, or should I be able to solve the real equation with Mathematica?
Thank you.
Edit: This is getting weirder. I ran A=1; tau = NDSolve[{A*Th''[r] == Th[r]^4, Th'[1] == 0, Th[0] == 1}, Th, {r, 0, 1}]; T[x_] = Evaluate[Th[x]/.tau]; Plot[T[x], {x, 0, 1}] for different values of A, and NDSolve only crashes for values smaller than A = 0.4821. I have no idea where to go from there. The correct value for A is A = 0.001356.
- As you receive help, try to give it too, by answering questions in your area of expertise.
- Take the tour and check the faqs!
- When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign!
– Mar 09 '16 at 15:34you will see what the solution looks like. $T'(0) < 0$ in this solution, which is correct. However, the real values do not match the physical units in this reduced solution.
– RegencyAndCo Mar 09 '16 at 20:49