Documentation |
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".
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[\mathrm{0..1}\right]\end{array}$$
$$Dg\left(0\right)=g\left(1\right)=0$$
Seek an approximate solution by collocation from C^{1}^{ } 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;
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(:).';
With this, the numerical problem you want to solve is to find $$y\in {S}_{4,knots}$$ that satisfies the nonlinear system
$$\begin{array}{c}Dy(0)=0\\ {(}^{y}+\epsilon {D}^{2}y(x)=1\text{for}x\text{}\in \text{colsites}\\ y(1)=0\end{array}$$
If y is your current approximation to the solution, then the linear problem for the supposedly better solution z by Newton's method reads
$$\begin{array}{c}Dz(0)=0\\ {w}_{0}(x)z(x)+\epsilon {D}^{2}z(x)=b(x)\text{for}x\text{}\in \text{colsites}\\ z\text{(1)=0}\end{array}$$
with w_{0}(x)=2y(x),b(x)=(y(x))^{2}+1. In fact, by choosing
$$\begin{array}{l}{w}_{0}(1):=1,\text{}{w}_{1}(0):=1\\ {w}_{1}(x):=0,\text{}{w}_{2}(x):=\epsilon \text{for}x\in \text{colsites}\end{array}$$
and choosing all other values of w_{0},w_{1},w_{2}, b not yet specified to be zero, you can give your system the uniform shape
$$\begin{array}{ccc}{w}_{0}\left(x\right)z\left(x\right)+{w}_{1}\left(x\right)Dz\left(x\right)+{w}_{2}\left(x\right){D}^{2}z\left(x\right)=b\left(x\right),& \text{for}& x\text{}\in \text{sites}\end{array}$$
with
sites = [0,colsites,1];
Because z∊S_{4,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 x∊sites 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 w_{0}(x),w_{1}(x),w_{2}(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 x^{2}–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
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 z–y 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.