6
$\begingroup$

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.

Improve this question
$\endgroup$
3
  • $\begingroup$ Welcome! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) 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! $\endgroup$
    – user9660
    Commented Mar 9, 2016 at 15:34
  • $\begingroup$ Do you have any a priori idea what the solution should look like? Since it's a boundary value problem, you might look into Mathematica's Shooting method: maybe you could give better initial guesses for Th'[0]. $\endgroup$
    – Chris K
    Commented Mar 9, 2016 at 20:28
  • $\begingroup$ If you run tau = NDSolve[{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}, PlotRange -> All] you 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. $\endgroup$
    – RegencyAndCo
    Commented Mar 9, 2016 at 20:49

2 Answers 2

Reset to default
4
$\begingroup$

When using Method->"Shooting" in NDSolve it helps to give good initial guesses. For example,

A = 0.001356;
tau = NDSolve[{A*Th''[r] == Th[r]^4, Th'[1] == 0, Th[0] == 1}, Th, {r, 0, 1}, 
Method -> {"Shooting", "StartingInitialConditions" -> {Th[0] == 1, Th'[0] == -17.1747}}];

seems to work fine.

Your problem seems very sensitive to that initial slope. How'd I get that crazy guess? I started with a solvable A value and then decreased it, using linear extrapolation to update the initial guess. Since A seems to vary over orders of magnitude, I stepped through it in fractional powers of 10.

Clear[A];
DTh0 = -0.5; (* initial guess *)
DTh = -0.5; (* "" *)
dAp = -0.01; (* A power step size *)
res = {};
Do[
  A = 10.^Ap;
  tau = NDSolve[{A*Th''[r] == Th[r]^4, Th'[1] == 0, Th[0] == 1}, Th, {r, 0, 1}, 
  Method -> {"Shooting", "StartingInitialConditions" -> {Th[0] == 1, Th'[0] == DTh0}}][[1]];

  DThold = DTh;
  DTh = Th'[0] /. tau;
  DTh0 = 2 DTh - DThold; (* linear extrapolation for next guess *)

  AppendTo[res, {A, DTh}];
,{Ap, 0, -3, dAp}]

Then to get the guess for your A value,

Interpolation[res][0.001356]
(* -17.1747 *)

It's sort of a hack of a solution. Hopefully someone can offer a more sophisticated approach!

Improve this answer
$\endgroup$
2
  • $\begingroup$ This is amazing, thank you. Now my answer is off, because this is a heat in = heat out problem such that $-kT'(0) = \int_0^L \epsilon \sigma T^4$, which is not true with the solution I just got, but now I can actually focus on the physics to sort this out! $\endgroup$
    – RegencyAndCo
    Commented Mar 10, 2016 at 7:09
  • $\begingroup$ Cannot edit my previous comment, but the solution is in fact correct, because the balance is actually $-ktT'(0) = \int_0^L \epsilon \sigma T^4$. All good. $\endgroup$
    – RegencyAndCo
    Commented Mar 10, 2016 at 9:11
4
$\begingroup$

Because the ODE is autonomous (i.e., r does not appear explicitly in it), it is "almost" possible to solve it symbolically, starting with DSolve with boundary conditions omitted. (DSolve returns unevaluated, if boundary conditions are included.)

DSolve[A*Th''[r] == Th[r]^4, Th, r]
(* Solve[(Hypergeometric2F1[1/5, 1/2, 6/5, -((2 Th[r]^5)/(5 A C[1]))]^2 Th[
r]^2 
   (5 C[1] + (2 Th[r]^5)/A))/(C[1]^2 (5 + (2 Th[r]^5)/(A C[1]))) == (r + C[2])^2, Th[r]] *)

To determine the two constants, extract and simply the resulting equation.

Simplify[First@%]
(* (Hypergeometric2F1[1/5, 1/2, 6/5, -((2 Th[r]^5)/(5 A C[1]))]^2 Th[r]^2)/C[1] == 
   (r + C[2])^2 *)

Next, take the square root of both sides of the equation, for which a positive or negative sign must be assumed for r + C[2]. Choose the latter.

Reverse[Simplify[Sqrt[#], r + C[2] < 0 && 
    Hypergeometric2F1[1/5, 1/2, 6/5, -((2 Th[r]^5)/(5 A C[1]))] Th[r] > 0] & /@ %]
(* -r - C[2] == Sqrt[1/C[1]]
   Hypergeometric2F1[1/5, 1/2, 6/5, -((2 Th[r]^5)/(5 A C[1]))] Th[r] *)

Applying the Th[0] == 1 boundary condition readily determines C[2].

Rule @@ (% /. r -> 0 /. Th[0] -> 1);
eq1 = (%% /. %) /. C[1] -> -2 c/(5 A)
(* -r + Sqrt[5/2] Sqrt[-(A/c)] Hypergeometric2F1[1/5, 1/2, 6/5, 1/c] == 
   Sqrt[5/2] Sqrt[-(A/c)]Hypergeometric2F1[1/5, 1/2, 6/5, Th[r]^5/c] Th[r] *)

where the substitution C[1] -> -2 c/(5 A) has been introduced for simplicity. Next, apply the Th'[1] == 0 boundary condition.

FullSimplify[D[eq1, r]] /. r -> 1
(* 2 + (Sqrt[10] Sqrt[-(A/c)] Th'[1])/Sqrt[1 - Th[1]^5/c] == 0 *)

which cannot be satisfied unless c == Th[1]^5.

eq2 = eq1 /. c -> Th[1]^5
(* -r + Sqrt[5/2]Hypergeometric2F1[1/5, 1/2, 6/5, 1/Th[1]^5] Sqrt[-(A/Th[1]^5)] == 
   Sqrt[5/2] Hypergeometric2F1[1/5, 1/2, 6/5, Th[r]^5/Th[1]^5] Sqrt[-(A/Th[1]^5)] Th[r] *)

eq2 plus the boundary condition eq2/.r -> 0 uniquely determine Th[r]. However, Th[1] must be computed numerically, here for A == 1356 10^-6.

FindRoot[eq2 /. r -> 1 /. A -> 1356 10^-6, {Th[1], 1/10}, MaxIterations -> 1000] // Chop
eq3 = eq2 /. % /. A -> 1356 10^-6
(* {Th[1] -> 0.137621} *)
(* (1. + 1.42981 I) - r == 
   (0. + 8.28684 I) Hypergeometric2F1[1/5, 1/2, 6/5, 20257.2 Th[r]^5] Th[r] *)

which can be plotted parametrically. (Note that, if r + C[2] > 0 had been assumed when square roots were taken above, eq2 now would have no real zeroes for Th[1]).

temin = Th[1] /. %%;
arg = {Subtract @@ eq3 /. Th[r] -> tem /. r -> 0, tem};
ParametricPlot[arg // Chop, {tem, temin, 1}, AspectRatio -> 1/GoldenRatio, 
    PlotRange -> {{0, 1}, {.1, 1}}, AxesLabel -> {r, Th}, 
    LabelStyle -> Directive[Black, Bold, 14], ImageSize -> Large]

enter image description here

Th[1] as a function of A can be plotted in a similar way.

Equal @@ (Solve[eq2 /. r -> 1, A][[1, 1]])
(* A == -((2 Gamma[7/10]^2 Th[1]^5)/(5 (Gamma[7/10] 
   Hypergeometric2F1[1/5, 1/2, 6/5, 1/Th[1]^5] - Sqrt[π] Gamma[6/5] Th[1])^2)) *)

arg = {Last[%] /. Th[1] -> tem, tem};
ParametricPlot[arg // Chop, {tem, 10^-3, 1/2}, AspectRatio -> 1/GoldenRatio, 
    AxesLabel -> {A, "Th[1]"}, LabelStyle -> Directive[Black, Bold, 14], 
    ImageSize -> Large, Epilog -> {Red, PointSize[Large], Point[{1356 10^-6, temin}]}]

enter image description here

Values of {Th[1], A} for the first plot are superimposed as a red dot. That it lies on the steep part of the curve may explain the difficulty encountered by the OP in solving the ODE numerically.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.