Finite Difference Method – Nonlinear ODE

Section 4.5 Finite Difference Method – Nonlinear ODE

Subsection 4.5.1 Heat conduction with radiation

If we again consider the heat in a metal bar of length \(L\text{,}\) but this time consider the effect of radiation as well as conduction, then the steady state equation has the form

\begin{equation*} u_{xx}- d (u^{4} - u_{b}^{4}) = - g(x)\text{,} \end{equation*}

where \(u_{b}\) is the temperature of the background, \(d\) incorporates a coefficient of radiation and \(g(x)\) is the heat source.

🔗

If we again replace the continuous problem by its discrete approximation then we get

\begin{equation} \frac{u_{i+1}- 2 u_{i} + u_{i-1}}{h^{2}}- d(u_{i}^{4} - u_{b}^{4}) = - g_{i} = -g(x_{i})\text{.}\tag{4.5.1} \end{equation}

This equation is nonlinear in the unknowns, thus we no longer have a system of linear equations to solve, but a system of nonlinear equations. One way to solve these equations would be by the multivariable Newton method in Section 2.6. Instead, we introduce another iterative method.

🔗

Subsection 4.5.2 Relaxation Method for Nonlinear Finite Differences

We can rewrite equation (4.5.1) as

\begin{equation*} u_{i+1}- 2 u_{i} + u_{i-1}= h^{2} d(u_{i}^{4} - u_{b}^{4}) - h^{2} g_{i}\text{.} \end{equation*}

From this we can solve for \(u_{i}\) in terms of the other quantities:

\begin{equation*} 2 u_{i} = u_{i+1}+ u_{i-1}- h^{2} d(u_{i}^{4} - u_{b}^{4}) + h^{2} g_{i}\text{.} \end{equation*}

Next we add \(u_{i}\) to both sides of the equation to obtain

\begin{equation*} 3 u_{i} = u_{i+1}+ u_{i} + u_{i-1}- h^{2} d(u_{i}^{4} - u_{b}^{4}) + h^{2} g_{i}\text{,} \end{equation*}

and then divide by 3 to get

\begin{equation*} u_{i} = \frac{1}{3}\left( u_{i+1}+ u_{i} + u_{i-1}\right) - \frac{h^{2}}{3}\left( d(u_{i}^{4} - u_{b}^{4}) - g_{i} \right)\text{.} \end{equation*}

🔗

Now for the main idea, we will begin with an initial guess for the value of \(u_{i}\) for each \(i\text{,}\) which we can represent as a vector \(\ub^{0}\text{.}\) Then we will use the above equation to get better estimates, \(\ub^{1}\text{,}\) \(\ub^{2}\text{,}\) …, and hope that they converge to the correct answer.

🔗

If we let

\begin{equation*} \ub^{j} = (u_{0}^{j}, u_{1}^{j}, u_{2}^{j}, \ldots, u_{n-1}^{j}, u_{n}^{j}) \end{equation*}

denote the \(j\)th approximation, then we can obtain that \(j+1\)st estimate from the formula

\begin{equation*} u^{j+1}_{i} = \frac{1}{3}\left( u^{j}_{i+1}+ u^{j}_{i} + u^{j}_{i-1}\right) - \frac{h^{2}}{3}\left( d((u^{j}_{i})^{4} - u_{b}^{4}) - g_{i} \right)\text{.} \end{equation*}

Notice that \(g_{i}\) and \(u_{b}\) do not change. In the resulting equation, we have \(u_{i}\) at each successive step depending on its previous value and the equation itself.

🔗

Subsection 4.5.3 Implementing the Relaxation Method

In the following program we solve the finite difference equations ((4.5.1)) with the boundary conditions \(u(0) = 0\) and \(u(L) = 0\text{.}\) We let \(L = 4\text{,}\) \(n=4\text{,}\) \(d = 1\text{,}\) and \(g(x) = \sin(\pi x/4)\text{.}\) Notice that the vector u always contains the current estimate of the values of \(\ub\text{.}\)

🔗

Program 4.5.1. mynonlinheat

🔗

% mynonlinheat (lacks comments)
% Purpose:
L = 4;             %
n = 4;             %
h = L/n;           %
hh = h^2/3;        %
u0 = 0;            %
uL = 0;            %
ub = .5;           % 
ub4 = ub^4;        %
x = 0:h:L;         %
g = sin(pi*x/4);   %
u = zeros(1,n+1);  %
steps = 4;         %
u(1)=u0;           % 
u(n+1)=uL;         %
for j = 1:steps
     %
     u(2:n) = (u(3:n+1)+u(2:n)+u(1:n-1))/3 + hh*(-u(2:n).^4+ub4+g(2:n));
end
plot(x,u)

If you run this program with the given n and steps the result will not seem reasonable.

🔗

We can plot the initial guess by adding the command plot(x,u) right before the for loop. We can also plot successive iterations by moving the last plot(x,u) before the end. Now we can experiment and see if the iteration is converging. Try various values of \(steps\) and \(n\) to produce a good plot. You will notice that this method converges quite slowly. In particular, as we increase \(n\text{,}\) we need to increase \(steps\) like \(n^{2}\text{,}\) i.e. if \(n\) is large then steps needs to be really large.

🔗

Exercises 4.5.4 Exercises

1.

Modify the script program mynonlinheat (Program 4.5.1) to plot the initial guess and all intermediate approximations. Add complete comments to the program. Print the program and a plot using \(n = 12\) and steps large enough to see convergence.

🔗

2.

Modify your improved mynonlinheat to mynonlinheattwo that has the boundary conditions

\begin{equation*} u(0) = 5 \quad\text{and}\quad u(L) = 10\text{.} \end{equation*}

Fix the comments to reflect the new boundary conditions. Print the program and a plot using \(n= 50\) and large enough steps to see convergence.

🔗

Prev Top Next