## Solving A Nonlinear ODE

This section discusses these aspects of a nonlinear ODE problem:

You can run this example: “Solving a Nonlinear ODE with a Boundary Layer by Collocation”.

### Problem

Consider the nonlinear singularly perturbed problem:

`$\begin{array}{ccc}\epsilon {D}^{2}g\left(x\right)+{\left(g\left(x\right)\right)}^{2}=1& on& \left[0..1\right]\end{array}$`
`$Dg\left(0\right)=g\left(1\right)=0$`

### Approximation Space

Seek an approximate solution by collocation from C1 piecewise cubics with a suitable break sequence; for instance,

```breaks = (0:4)/4; ```

Because cubics are of order 4, you have

```k = 4; ```

Obtain the corresponding knot sequence as

```knots = augknt(breaks,k,2); ```

This gives a quadruple knot at both 0 and 1, which is consistent with the fact that you have cubics, i.e., have order 4.

This implies that you have

```n = length(knots)-k; n = 10; ```

i.e., 10 degrees of freedom.

### Discretization

You collocate at two sites per polynomial piece, i.e., at eight sites altogether. This, together with the two side conditions, gives us 10 conditions, which matches the 10 degrees of freedom.

Choose the two Gaussian sites for each interval. For the standard interval [–0.5,0.5] of length 1, these are the two sites

`gauss = .5773502692*[-1/2; 1/2];`

From this, you obtain the whole collection of collocation sites by

```ninterv = length(breaks)-1; temp = ((breaks(2:ninterv+1)+breaks(1:ninterv))/2); temp = temp([1 1],:) + gauss*diff(breaks); colsites = temp(:).'; ```

### Numerical Problem

With this, the numerical problem you want to solve is to find $y\in {S}_{4,knots}$ that satisfies the nonlinear system

### Linearization

If y is your current approximation to the solution, then the linear problem for the supposedly better solution z by Newton's method reads

with w0(x)=2y(x),b(x)=(y(x))2+1. In fact, by choosing

and choosing all other values of w0,w1,w2, b not yet specified to be zero, you can give your system the uniform shape

with

```sites = [0,colsites,1]; ```

### Linear System to Be Solved

Because zS4,knots, convert this last system into a system for the B-spline coefficients of z. This requires the values, first, and second derivatives at every xsites and for all the relevant B-splines. The command `spcol` was expressly written for this purpose.

Use `spcol` to supply the matrix

```colmat = ... spcol(knots,k,brk2knt(sites,3)); ```

From this, you get the collocation matrix by combining the row triple of `colmat` for x using the weights w0(x),w1(x),w2(x) to get the row for x of the actual matrix. For this, you need a current approximation y. Initially, you get it by interpolating some reasonable initial guess from your piecewise-polynomial space at the `sites`. Use the parabola x2–1, which satisfies the end conditions as the initial guess, and pick the matrix from the full matrix `colmat`. Here it is, in several cautious steps:

```intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:); coefs = intmat\[0 colsites.*colsites-1 0].'; y = spmak(knots,coefs.');```

Plot the initial guess, and turn hold on for subsequent plotting:

```fnplt(y,'g'); legend('Initial Guess (x^2-1)','location','NW'); axis([-0.01 1.01 -1.01 0.01]); hold on```

### Iteration

You can now complete the construction and solution of the linear system for the improved approximate solution z from your current guess y. In fact, with the initial guess y available, you now set up an iteration, to be terminated when the change zy is small enough. Choose a relatively mild ε = .1.

```tolerance = 6.e-9; epsilon = .1; while 1 vtau = fnval(y,colsites); weights=[0 1 0; [2*vtau.' zeros(n-2,1) repmat(epsilon,n-2,1)]; 1 0 0]; colloc = zeros(n,n); for j=1:n colloc(j,:) = weights(j,:)*colmat(3*(j-1)+(1:3),:); end coefs = colloc\[0 vtau.*vtau+1 0].'; z = spmak(knots,coefs.'); fnplt(z,'k'); maxdif = max(max(abs(z.coefs-y.coefs))); fprintf('maxdif = %g\n',maxdif) if (maxdif<tolerance), break, end % now reiterate y = z; end legend({'Initial Guess (x^2-1)' 'Iterates'},'location','NW');```

The resulting printout of the errors is:

```maxdif = 0.206695 maxdif = 0.01207 maxdif = 3.95151e-005 maxdif = 4.43216e-010```

If you now decrease ε, you create more of a boundary layer near the right endpoint, and this calls for a nonuniform mesh.

Use `newknt` to construct an appropriate finer mesh from the current approximation:

```knots = newknt(z, ninterv+1); breaks = knt2brk(knots); knots = augknt(breaks,4,2); n = length(knots)-k; ```

From the new break sequence, you generate the new collocation site sequence:

```ninterv = length(breaks)-1; temp = ((breaks(2:ninterv+1)+breaks(1:ninterv))/2); temp = temp([1 1], :) + gauss*diff(breaks); colpnts = temp(:).'; sites = [0,colpnts,1]; ```

Use `spcol` to supply the matrix

```colmat = spcol(knots,k,sort([sites sites sites])); ```

and use your current approximate solution `z` as the initial guess:

```intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:); y = spmak(knots,[0 fnval(z,colpnts) 0]/intmat.'); ```

Thus set up, divide ε by 3 and repeat the earlier calculation, starting with the statements

```tolerance=1.e-9; while 1 vtau=fnval(y,colpnts); . . . ```

Repeated passes through this process generate a sequence of solutions, for ε = 1/10, 1/30, 1/90, 1/270, 1/810. The resulting solutions, ever flatter at 0 and ever steeper at 1, are shown in the example plot. The plot also shows the final break sequence, as a sequence of vertical bars. To view the plots, run the example “Solving a Nonlinear ODE with a Boundary Layer by Collocation”.

In this example, at least, `newknt` has performed satisfactorily.