# Shooting Method with RK4 and Thermal Radiation

I am attempting to numerically solve the following problem. I decompose it into a system of two first order ODEs and then solve via the shooting method. I use the fourth order Runge-Kutta (RK4) method to solve each iteration of the shooting method.

$$k \frac{d^2T}{dx^2}=h(T-T_{\inf}) + \sigma(T^4-T^4_{surr})\\ T_{inf} = T_{surr} = 200 \, K\\ k = 1 \quad h = 0.05 \, \text{m}^{-2}\quad \sigma = 2.7*10^{-9} \, \text{K}^{-3}\,\text{m}^{-2}\quad T_L = 300 \, K \quad T_R = 400 \, K\\ L = 10 \, \text{m}$$

If I set $$\sigma = 0$$ then the solution converges for RK4 in just 3 iterations. I have checked it against the analytical solution for when $$\sigma=0$$.

However, if I keep the radiation term ($$\sigma\neq0$$) then RK4 quickly diverges during the first iteration of the shooting method. About half-way thru the solution domain it explodes.

Is there something special or unique about the following system that I am not aware? Some special case or sensitivity to initial conditions?

The system of first order ODEs that I am applying RK4 to is: $$\frac{dz}{dx} = w\\ \frac{dw}{dx} = h(z-T_{inf}) + \sigma (z^4-T^4_{inf})\\ \text{where} \quad z(0) = T_L\\ \text{and} \quad w(0) = 1 \; \text{(Guess 1)}$$

• Shooting isn't stable. I would recommend using a MIRK-based method instead in order to get a more stable rootfinding problem. The tutorials of DifferentialEquations.jl show how to quickly swap out a shooting problem for a MIRK problem. Feb 8 '20 at 21:26
• Thank you, I'll check out those resources, but this is an assignment for a class and the instructions are "use shooting method with RK4" Feb 9 '20 at 1:17
• Can you give the numerical values of all the parameters involved $h,\sigma,T_{inf},T_L$ and domain length. Feb 9 '20 at 16:48
• In the second system, at some point you have to divide by $k$, you can decide on the first or second equation. Is this only missing here or also in your code? Feb 9 '20 at 17:37
• It's not very effective to use the shooting method for a two-point boundary value problem. Just solve it as a coupled problem in space, using the finite difference or finite element method! Feb 10 '20 at 21:57