NONLINEAR BVP The temperature distribution along a fin can be modeled as a BVP [latex]T_{x x}-\beta T^{4}-\alpha T+\gamma=0, \quad T(0)=500, \quad T(0.2)=350, \quad T_{x x}=\frac{d^{2} T}{d x^{2}}[/latex] where [latex]\alpha=20, \beta=10^{-8}[/latex], and [latex]\gamma=5 \times 10^{3}[/latex] with all parameters in consistent physical units. Solve the nonlinear BVP using the shooting method.

There are two state variables selected [latex]\eta_{1}=T, \eta_{2}=T_{x}[/latex], and the state-variable equations are formed as [latex]\begin{cases}\eta_{1}^{\prime}=\eta_{2} & \eta_{1}(0)=500 \\ \eta_{2}^{\prime}=\alpha \eta_{1}+\beta \eta_{1}^{4}-\gamma^{\prime}, & \eta_{2}(0)=?\end{cases}[/latex] (8.3) Of the two required initial conditions, only [latex]\eta_{1}(0)=500[/latex] is available. As a first guess, we arbitrarily choose [latex]\eta_{2}(0)=100[/latex] so that Equation 8.3 becomes [latex] \begin{cases}\eta_{1}^{\prime}=\eta_{2} & \eta_{1}(0)=500 \\ \eta_{2}^{\prime}=\alpha \eta_{1}+\beta \eta_{1}^{4}-\gamma^{\prime} & \eta_{2}(0)=100\end{cases} [/latex] We will solve this system via RK4 using 100 points.

Question 8.2: NONLINEAR BVP The temperature distribution along a fin can b...

Numerical Methods for Engineers and Scientists Using MATLAB®

NONLINEAR BVP

The temperature distribution along a fin can be modeled as a BVP

$T_{x x}-\beta T^{4}-\alpha T+\gamma=0, \quad T(0)=500, \quad T(0.2)=350, \quad T_{x x}=\frac{d^{2} T}{d x^{2}}$

where $\alpha=20, \beta=10^{-8}$ , and $\gamma=5 \times 10^{3}$ with all parameters in consistent physical units. Solve the nonlinear BVP using the shooting method.

The 'Blue Check Mark' means that either the MATLAB code/script/answer provided in the answer section has been tested by our team of experts; or the answer in general has be fact checked.

Learn more on how do we answer questions.

There are two state variables selected $\eta_{1}=T, \eta_{2}=T_{x}$ , and the state-variable equations are formed as

$\begin{cases}\eta_{1}^{\prime}=\eta_{2} & \eta_{1}(0)=500 \\ \eta_{2}^{\prime}=\alpha \eta_{1}+\beta \eta_{1}^{4}-\gamma^{\prime}, & \eta_{2}(0)=?\end{cases}$ (8.3)

Of the two required initial conditions, only $\eta_{1}(0)=500$ is available. As a first guess, we arbitrarily choose $\eta_{2}(0)=100$ so that Equation 8.3 becomes

$\begin{cases}\eta_{1}^{\prime}=\eta_{2} & \eta_{1}(0)=500 \\ \eta_{2}^{\prime}=\alpha \eta_{1}+\beta \eta_{1}^{4}-\gamma^{\prime} & \eta_{2}(0)=100\end{cases}$

We will solve this system via RK4 using 100 points.

>> f = @(x,eta)([eta(2);20*eta(1)+1e-8*eta(1)^4-5e3]);
>> x = linspace(0,0.2);
>> eta0_first = [500;100];
>> eta_first = RK4System(f,x,eta0_first);
>> T_end_first = eta_first(1,end)

T_end_first =

646.2081

The result overshoots the target of 350. Therefore, we must pick a second guess for $\eta_{2}(0)$ that leads to a value below the target. This will require at least one trial. It turns out that a second guess of $\eta_{2}(0)=-2000$ will work here. Then, Equation 8.3 reduces to

$\begin{cases}\eta_{1}^{\prime}=\eta_{2} & \eta_{1}(0)=500 \\ \eta_{2}^{\prime}=\alpha \eta_{1}+\beta \eta_{1}^{4}-\gamma^{\prime} & \eta_{2}(0)=-2000\end{cases}$

Note that $f$ and $x$ in our MATLAB code remain unchanged; only the initial state vector is modified.

>> eta0_second = [500;-2000];
>> eta_second = RK4System(f,x,eta0_second);
>> T_end_second = eta_second(1,end)

T_end_second =

157.0736

The result is below the target of 350 . In summary,

$\begin{array}{ll} \eta_{2}(0)=100 & T(0.2)=646.2081 \\ \eta_{2}(0)=-2000 & T(0.2)=157.0736 \end{array}$

Therefore, the unique value of $\eta_{2}(0)$ that leads to the true solution lies in the interval [-2000, 100]. We will find this value using the bisection method. The iterations (maximum 30) terminate when the computed $T(0.2)$ is within $10^{-4}$ of the target.

eta20L = -2000; eta20R = 100; % Left and right end of the interval
kmax = 30; % Maximum number of iterations
tol = 1e-4; % Tolerance
T_end = 350; % Boundary condition (target)

for k = 1:kmax,
eta20 = (eta20L + eta20R)/2; % Bisection
eta0 = [500;eta20]; % Set initial state vector
eta = RK4System(f,x,eta0); % Solve the system
T(k) = eta(1,end); % Extract T(0.2)
err = T(k) - T_end; % Compare with target

% Adjust the left or right value of initial condition based on whether
% error is positive or negative

if abs(err) < tol,
break
end
if err > 0,
eta20R = eta20;
else
eta20L = eta20;
end
end

>> k

k =

21 % Number of iterations needed to meet tolerance

>> T(21)

ans =

350.0000 % Agrees with target T(0.2)=350

>> eta20

eta20 =

-1.1643e+03 % Actual value for the missing initial condition

With this information, Equation 8.3 becomes

$\begin{cases}\eta_{1}^{\prime}=\eta_{2} & \eta_{1}(0)=500 \\ \eta_{2}^{\prime}=\alpha \eta_{1}+\beta \eta_{1}^{4}-\gamma^{\prime}, & \eta_{2}(0)=-1164.3\end{cases}$

This system is then solved via RK4. The numerical solution $T$ of the original BVP is obtained by extracting the first row of the output, which represents the first state variable $\eta_{1}=T$ . Similarly, the first row of eta_first gives the solution corresponding to the first guess of $\eta_{2}(0)$ , and the first row of eta_second gives the solution corresponding to the second guess of $\eta_{2}(0)$ . This is summarized in the script below.

>> T_first = eta_first(1,:);
>> T_second = eta_second(1,:);
>> eta0 = [500;eta20];
>> eta_final = RK4System(f,x,eta0);
>> T_final = eta_final(1,:);
>> plot(x,T_first,x,T_second,x,T_final) % Figure 8.3

The shooting method loses its efficiency when applied to higher-order boundaryvalue problems, which will require more than one guess for the initial values. For those cases, other techniques, such as the finite-difference method, need to be employed.