Question 6.4: Consider the same ODE as in Example 6.3 but use a Neumann bo...

Numerical Methods with Chemical Engineering Applications

Consider the same ODE as in Example 6.3 but use a Neumann boundary condition on the left boundary, e.g.,

$\frac{dy}{dx}|_{x=−1}= 10$ (6.3.29)

In this case, the form of the equations for the interior nodes and the right boundary are the same as before, given by Eqs. (6.3.16) and (6.3.18). For the boundary condition with the fictitious node, we get
$\frac{y_2 − y_0}{2 Δx }= 10$ (6.3.30)

$\frac{y_{i+1} − 2y_i + y_{i−1}}{Δx^2} +2y_i(1 − 2y^2_i ) = 0$ (6.3.16)

$y_n = 1$ (6.3.18)

which we solve to get
$y_0 = y_2 - 20 Δx$ (6.3.31)
We now discretize the differential equation at $x_1$ , to get

$\frac{y_2 − 2y_1 + y_0}{Δx ^2} + 2y_1(1 − 2y^2_1) = 0$ (6.3.32)

Substituting for $y_0$ gives a new first row for the residual

$R_1 = \frac{2y_2 − 2y_1 − 20 Δx}{Δx ^2}+ 2y_1(1 − 2y^2_1) = 0$ (6.3.33)
We thus modify the program getR from the last example:

1 function R = getR(y,dx,n)
2 R = zeros(n,1);
3 R(1,1) = (2*y(2)-2*y(1)-20*dx)/dx^2 + 2*y(1)*(1-2*y(1)^2);
4 for i = 2:n-1
5 R(i,1) = (y(i+1)-2*y(i)+y(i-1))/(dx^2)+2*y(i)*(1-2*y(i)^2);
6 end
7 R(n,1) = y(n) - 1;

The only change is in line 3, where we now write Eq. (6.3.33) in place of the Dirichlet condition.
Since we have a new $R_1$ , the non-zero entries in the first row of the Jacobian are now

$J_{1,1} =\frac {−2}{Δx^2} + 2 − 12y^2_1$ (6.3.34)

$J_{1,2} = \frac{2}{Δx^2}$ (6.3.35)

We thus need to make small changes to the program to compute the Jacobian as well:

1 function J=getJ(y,dx,n)
2 J = zeros(n);
3 denom = dx^(-2);
4 J(1,1) = -2*denom + 2 - 12*y(1)^2;
5 J(1,2) = 2*denom;
6 for i = 2:n-1
7 J(i,i-1) = denom;
8 J(i,i) = -2*denom + 2 - 12*y(i)^2;
9 J(i,i+1) = denom;
10 end
11 J(n,n) = 1;

Lines 4 and 5 of this program compute the new entries in the Jacobian for the Neumann boundary condition.
These changes to the programs for the residual and Jacobian are all we need; the rest of the solution is identical to Example 6.3. If we also continue to use the same initial guess, it takes eight iterations for Newton–Raphson to converge to the solution in Fig. 6.6. The solution y(x) is now much more nonlinear due to the very steep boundary condition at x = −1.