Question 6.3: Develop a MATLAB code to solve d^2y/dx^2 = −2y(1 − 2y^2), y(...

Following the procedure above, for the interior nodes we have the following algebraic equations

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

For the boundary nodes, in this case of Dirichlet boundary conditions, we can directly fix the values of the variables y_1  and  y_n,

y_1 = −1            (6.3.17)
and
y_n = 1         (6.3.18)
The above can now be put in a residual form R(y) = 0 with

R_1 = y_1 + 1

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

R_n = y_n − 1                     (6.3.19)

for i = 2, ... , n − 1. To solve this problem with Newton–Raphson, we need to be able to generate the residual vector for an arbitrary number of nodes:

This function takes in the values of the unknowns, the spacing between nodes, and the number of nodes. While the latter two variables are not independent, we chose to send them both in this example rather than compute, for example, Δx from n. On line 2, we make space in the memory for the residual. If we use many nodes, we can get a major slowdown in MATLAB if we forget to initialize matrices and vectors to the correct size before filling them in! (In other structured programming languages, you generally have to declare all of your variables at the outset anyway, so this is a problem that really only arises when you program in MATLAB.) After making space in the memory, we fill in the vector with the residual equations. Line 3 handles the left boundary condition, lines 4–6 loop over the interior of the domain, and line 7 writes the entry for the right boundary.

We also need to specify the non-zero entries in the Jacobian matrix. For the first line, there is only one non-zero entry,
J_{1,1} = 1         (6.3.20)
Likewise, for the last row, there is also only one non-zero entry,
J_{n,n} = 1                  (6.3.21)

Note that we could reduce the problem from an n × n system to an (n − 2) × (n − 2) system by substituting for the first and last nodes. This is a matter of preference, since for n >> 1 the time saved by this substitution is not really significant. We will generally include the boundary conditions in all of the problems that we solve here, since the code is easier to relate back to the original differential equation. For the interior rows of the Jacobian, corresponding to the index i = 2, 3, ... , n−1, there are three non-zero entries

J_{i,i−1} =\frac{1}{Δx^2}                 (6.3.22)

J_{i,i} = \frac{−2}{Δx^2} + 2 − 12y^2_i                                            (6.3.23)

J_{i,i+1} =\frac{ 1}{Δx^2}                                  (6.3.24)

The MATLAB program to write the Jacobian is then:

The structure of the program is quite similar to getR; the function takes in the vector of unknowns and the discretization parameters, makes space in the memory to write the Jacobian, and then fills in the non-zero entries. Notice in line 4 that we have defined a new variable denom to be equal to ( Δx)^{−2}.  This step removes quite a lot of computation, since we calculate this term 3(n − 2) times when we write the Jacobian. There are often many opportunities to speed up your codes, but we generally suggest that you write the simplest code possible as you are learning the methods and worry about speed later. The reason for this advice is that most increases in speed will reduce the readability of the code and make it harder to find mistakes.

To solve this problem, we need to make an initial guess and then iterate using Newton–Raphson. Let’s go back to Program 3.5 and make a few small changes:

Obviously, we need to send the values of Δx and n to this program, since they are needed by getR and getJ. We again are using a banded solver due to the tridiagonal structure of the Jacobian. Note that the banded solver can be even more important for nonlinear problems than linear problems, since we may need to iterate a number of times before we converge on a solution to a nonlinear problem.
The only remaining item is to decide on an initial guess. Based on the boundary conditions, the choice is y^{(0)}_i = x_i since it satisfies the boundary conditions and the homogeneous problem for Eq. (6.3.15). As we see in Fig. 6.4, this guess is adequate to get the solution to converge. Indeed, the nonlinearity is relatively weak and the end result is quite similar to a straight line with some small fluctuations.