Solving one dimension steady state heat equation with finite differences

[Arkady87]

I have a project where I need to solve

T''(x) = bT^4 ; 0<=x<=1
T(0) = 1
T'(1) = 0

using finite differences to generate a system of equations in Matlab and solve the system to find the solution

So far I have:
(using centred 2nd degree finite difference)
T''(x) = (T(x+h) - 2T(x) + T(x-h)) / h^2 = bT(x)^4
and
(using 2 order backward difference)
T'(x) = (3T(x) - 4T(x-h) + T(x-2h)) / 2h

I just don't know how to write the code that will make the system of n non-linear equations to solve. I have a function that will solve them if I can get there...

any help or nudge in the right direction is appreciated

 

[haruspex]

Are you sure the equation is right? Doesn't look like any heat equation I know. It looks like a diffusion equation on the left and a radiation equation on the right.

 

[Chestermiller]

I have a project where I need to solve

T''(x) = bT^4 ; 0<=x<=1
T(0) = 1
T'(1) = 0

using finite differences to generate a system of equations in Matlab and solve the system to find the solution

So far I have:
(using centred 2nd degree finite difference)
T''(x) = (T(x+h) - 2T(x) + T(x-h)) / h^2 = bT(x)^4
and
(using 2 order backward difference)
T'(x) = (3T(x) - 4T(x-h) + T(x-2h)) / 2h

I just don't know how to write the code that will make the system of n non-linear equations to solve. I have a function that will solve them if I can get there...

any help or nudge in the right direction is appreciated

This is a steady state heat conduction/radiation problem. The non-linearity is a difficulty. You might try to solve it as a transient heat conduction/radiation problem by addition a time derivative to the equation in the appropriate place, and solving for how the temperature varies with time as x. Here, you would use the so-called Method of Lines. You might try an explicit finite differencing approach, and, with small enough time steps should get convergence. Alternately, you could try an implicit method (temperatures at forward time step), and use an automatic package ODE solver to do the integration. For this problem, the non-stiff option should be adequate. Try as the initial condition T = (1-x)^2

 

[haruspex]

Does this help?
T'' T' = b T' T⁴
(T')² = 2b T⁵/5 + c
x = ∫(2b T⁵/5 + c)^-1/2dT

 

[jackmell]

I just don't know how to write the code that will make the system of n non-linear equations to solve.

Yeah me neither but I got an idea before I try to find it in a book. How about you solve for both the values of [itex]t_i[/itex] and [itex]t_i^'[/itex] all the way down the line and for starters, just break it up into 10 intervals.

So you have:

[tex]t_{i+1}-2t_i+\frac{t_{i-1}}{0.1^2}=bt_i^4[/tex]

as well as compute the derivative at each point:

[tex]d_i=\frac{3t_i-4 t_{i-1}+t_{i-2}}{0.2}[/tex]

and then solve for the set [itex](t_i, d_i)[/itex]

Not sure though ok but it's a start and in math, starts are important. :)

 

[haruspex]

Assuming you just want to solve the equation numerically, whether using Matlab or otherwise...
x = ∫(2b t⁵/5 + c)^-1/2dt
Writing B² = 2b/5, T'(1) = 0 implies c = -B²α⁵, where α = T(1). T(0) = 1 gives:
[tex]x = B^{-1}∫^1_{t=T}(t^5 - α^5)^{-1/2}dt[/tex]
It remains to find [itex]\hat{α}[/itex] s.t.
[tex]B = ∫^1_{t=\hat{α}}(t^5 - \hat{α}^5)^{-1/2}dt[/tex]
Consider
[tex]I(α) = ∫^1_{t=α}(t^5 - α^5)^{-1/2}dt[/tex]
[tex]I(α/k) = ∫^1_{t=α/k}(t^5 - (α/k)^5)^{-1/2}dt = ∫^k_{t=α}((t/k)^5 - (α/k)^5)^{-1/2}d(t/k)[/tex]
[tex] = k^{3/2}∫^k_{t=α}(t^5 - α^5)^{-1/2}dt[/tex]
So we can proceed as follows:
Choose an arbitrary [itex]α, \check{α}[/itex] = 0.1 say.
Perform the integral [tex]F(k) = ∫^k_{t=\check{α}}(t^5 - \check{α}^5)^{-1/2}dt[/tex] stepwise on k. At each step, check whether F(k)k^3/2 exceeds B. When it does, backtrack one step and binary chop the last step size to refine the value of k. Then [itex]\hat{α} = \check{α}/k[/itex]
Note that starting the calculation of the integral is interesting because the integrand is infinite initially. This can be managed by approximating
[tex]((α+δt)^5 - α^5)^{-1/2} = (5α^4δt)^{-1/2}[/tex]
and integrating to produce
[tex]F(\check{α}+δt) = 2(5\check{α}^4δt)^{1/2}[/tex]