Question

! Required information An insulated heated rod with a uniform heat source can be modeled with the Poisson equation as given below. de = -f(x) Given: A heat source fx) = 25°C/m2 and the boundary conditions Tlx = 0) = 40 °C and Tlx = 10) = 200 °C Solve the ODE using the shooting method. (Round the final answer to the nearest integer.) The temperature distribution at x = 4 is K. Use the fourth-order RK method using a step size of 0.125.

Answer 1

Shooting method:

Suppose we have a second order ODE:

T" = f(2,1,T'), a< x <b

where a prime denotes a derivative with respect to x, with boundary conditions:

Ta) = a T(6) = 8

Then we consider an initial value problem:

T" = f(x,T,T'), a< x <b, T(a) = T(a) = a, T'(a) =

where s is a paramter given by:

$s = \frac{\beta - T_1(b)}{T_2(b)}$

where, is the solution of initial value problem:

T12

$T^{''} = f(x,T,T^{'})\,,\quad a<x<b\,,\quad\,T(a)=\alpha\,,\quad T^{'}(a)=0$

and $T_2(x)$ is the solution of initial value problem:

$T^{''} = f(x,T,T^{'})\,,\quad a<x<b\,,\quad\,T(a)=0\,,\quad T^{'}(a)=1$

Then, the solution of our original boundary value problem is:

Now, our boundary value problem (BVP) is:

$T^{''} = -f(x)\,,\quad 0<x<10\,,\quad T(0)=40\,,\quad T(10)=200$

so, $a = 0\,,\quad b = 10\,,\quad \alpha = 40\,,\quad \beta = 200$

We need to find the solutions to the initial value problems (IVP):

$T^{''} = -f(x)\,,\quad 0<x<10\,,\quad T(0)=40\,,\quad T^{'}(0)=0$

$T^{''} = -f(x)\,,\quad 0<x<10\,,\quad T(0)=0\,,\quad T^{'}(0)=1$

using fourth order Runge Kutta (RK4) method.

To use the RK4 method, we need to convert the second order ODE into two first order ODEs:

Let, , so that we have two first order ODEs:

U

$U^{'} = -f(x)\,,\quad 0<x<10\,,\quad U(0) = 0$

$T^{'} = U\,,\quad 0<x<10\,,\quad T(0) = 40$

Similarly, for the second IVP, we will have and

Now the RK4 scheme for solving a system of coupled ODEs:

$y'=f(x,y,z)\,,\quad y(0) = p$

$z'=g(x,y,z)\,,\quad z(0) = q$

is:

$k_1=h\,f(x_i+h/2,y_i+k_0/2,z_i+l_0/2)$

$l_1=h\,g(x_i+h/2,y_i+k_0/2,z_i+l_0/2)$

$k_2=h\,f(x_i+h/2,y_i+k_1/2,z_i+l_1/2)$

$l_2=h\,g(x_i+h/2,y_i+k_1/2,z_i+l_1/2)$

$k_3=h\,f(x_i+h,y_i+k_2,z_i+l_2)$

$l_3=h\,g(x_i+h,y_i+k_2,z_i+l_2)$

$y_{i+1} = y_i + \frac{1}{6}\,(k_0+2k_1+2k_2+k_3)$

$z_{i+1} = z_i + \frac{1}{6}\,(l_0+2l_1+2l_2+l_3)$

$y_0 = p\,,\quad z_0=q$

where, h is the step size in x

In this way, we solve the two IVPs and get s and then get the solution for our BVP.

The coding is done in python.

############## shooting_method.py #########

import numpy as np
import matplotlib.pyplot as plt

a = 0.0
b = 10.0
alpha = 40.0
beta = 200.0
f = 25.0
h = 0.125
x = np.arange(a,b+h,h)
nitr = len(x)
T1 = np.zeros(nitr)
T2 = np.zeros(nitr)
U1 = np.zeros(nitr)
U2 = np.zeros(nitr)
T1[0] = alpha
T2[0] = 0.0
U1[0] = 0.0
U2[0] = 1.0

for i in range(nitr-1):
# Solving for U1 and T1
k0 = h*(-f)
l0 = h*U1[i]
k1 = h*(-f)
l1 = h*(U1[i]+0.5*l0)
k2 = h*(-f)
l2 = h*(U1[i]+0.5*l1)
k3 = h*(-f)
l3 = h*(U1[i]+l2)
U1[i+1] = U1[i] + (k0 + 2.0*k1 + 2.0*k2 + k3)/6.0
T1[i+1] = T1[i] + (l0 + 2.0*l1 + 2.0*l2 + l3)/6.0

# Solving for U2 and T2
k0 = h*(-f)
l0 = h*U2[i]
k1 = h*(-f)
l1 = h*(U2[i]+0.5*l0)
k2 = h*(-f)
l2 = h*(U2[i]+0.5*l1)
k3 = h*(-f)
l3 = h*(U2[i]+l2)
U2[i+1] = U2[i] + (k0 + 2.0*k1 + 2.0*k2 + k3)/6.0
T2[i+1] = T2[i] + (l0 + 2.0*l1 + 2.0*l2 + l3)/6.0
# for loop ends here

s = (beta - T1[nitr-1])/T2[nitr-1]
T = T1 + s*T2
plt.plot(x,T)
plt.xlabel('x')
plt.ylabel('Temperature in degC')
plt.title('T vs x solved using shooting method')
plt.grid()
plt.savefig('shooting-T_vs_x.png')
plt.show()
idx = np.where(x == 4.0)[0][0]
ans = T[idx] # in degC
ansK = ans + 273.15
mystr = 'Tmperature distribution at x = 4 is '+str(round(ansK,2))+' K'
print(mystr)

################ end of code ###########

screenshots of the code:

import numpy as np import matplotlib.pyplot as plt a = 0.0 b = 10.0 alpha = 40.0 beta = 200.0 f = 25.0 h = 0.125 X = np.arange(a,b+h, h) nitr = len(x) T1 = np.zeros(nitr) T2 = np.zeros(nitr) ui np.zeros(nitr) U2 = np.zeros(nitr) Ti[0] = alpha T2[0] = 0.0 U1[0] = 0.0 U2[@] = 1.0

for i in range(nitr-1): # Solving for ul and ti ko = h*(-f) 10 = h*11[i] kl = h*(-f) 11 h* (U1[i]+0.5*10) k2 : h*(-f) 12 h*(U1[i]+0.5*11) k3 h*(-f) 13 = h* (U1[i]+12) 11[i+1] = U1[i] + (ko + 2.0*k1 + 2.0*k2 + k3)/6.0 T1[i+1] = T1[i] + (10 + 2.0*11 + 2.0*12 + 13)/6.0 # Solving for U2 and T2 ko = :h*(-f) 10 = h*U2[i] ki : h*(-f) li h* (U2[i]+0.5*10) k2 h*(-f) 12 h* (U2[i]+0.5*11) k3 = h*(-f) 13 h*(U2[i]+12) U2[i+1] : U2[i] + (ko + 2.0*k1 + 2.0*k2 + k3)/6.0 T2[i+1] T2[i] + (10 + 2.0*11 + 2.0*12 + 13)/6.0 # for loop ends here

S = (beta - Ti[nitr-1])/T2[nitr-1] T = 11 + 5*T2 plt.plot(x,T) plt.xlabel('x') plt.ylabel('Temperature in degc') plt.title('T vs x solved using shooting method') plt.grid() plt.savefig('shooting-T_vs_x.png') plt.show() idx = np.where(x 4.0)[0][0] ans = T[idx] # in degC ansk = ans + 273.15 mystr = 'Tmperature distribution at x = 4 is '+str(round(ansk,2))+' K' print(mystr)

plot of the solution:

T vs x solved using shooting method 200 180 160 140 Temperature in degc 120 100 80 60 40 4 8 10 х

output on screen:

$ python3 shooting method.py Tmperature distribution at x = 4 is 335.34 K $

So, upto the nearest integer, the temperature distributionat x = 4 is 335 K