A Nonlinear ODE :: Advanced Examples (Spline Toolbox™)

Learn more about Spline Toolbox

A Nonlinear ODE

This section discusses these aspects of a nonlinear ODE problem:

The example can be run via the demo "Solving a Nonlinear ODE with a Boundary Layer by Collocation".

Problem

We consider the nonlinear singularly perturbed problem:

Approximation Space

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

breaks = (0:4)/4;

Since cubics are of order 4, we have

k = 4;

We 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 we have cubics, i.e., have order 4.

This implies that we have

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

i.e., 10 degrees of freedom.

Discretization

We 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.

We choose the two Gaussian sites for each interval. For the standard interval [-.5,.5] of length 1, these are the two sites

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

From this, we 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 we want to solve is to find that satisfies the nonlinear system

Linearization

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

with , . In fact, by choosing

and choosing all other values of not yet specified to be zero, we can give our system the uniform shape

with

sites = [0,colsites,1];

Linear System to Be Solved

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

We use spcol to supply the matrix

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

From this, we get the collocation matrix by combining the row triple of colmat for using the weights to get the row for of the actual matrix. For this, we need a current approximation . Initially, we get it by interpolating some reasonable initial guess from our piecewise-polynomial space at the sites. We use the parabola , 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.');

Iteration

We can now complete the construction and solution of the linear system for the improved approximate solution z from our current guess y. In fact, with the initial guess y available, we now set up an iteration, to be terminated when the change is small enough. We choose a relatively mild .

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 we now decrease , we create more of a boundary layer near the right endpoint, and this calls for a nonuniform mesh.

We 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, we 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];

We use spcol to supply the matrix

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

and use our 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, we cut 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 plot above. The plot also shows the final break sequence, as a sequence of vertical bars.

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

Least-Squares Approximation by "Natural" Cubic Splines

Construction of the Chebyshev Spline

Recommended Products

Get MATLAB trial software

Includes the most popular MATLAB recorded presentations with Q&A sessions led by MATLAB experts.

Get the Interactive Kit