1. The original problem
The Carrier equation can be found in Section 9.7 of the well-known textbook of Bender & Orszag [1], and is given as follows:
$$ \varepsilon u'' + 2(1-x^2)u + u^2 = 1,\qquad u(-1) = 0,~~ u(1) = 0. $$
This is a nonlinear ODE boundary-value problem with multiple solutions. We can find a solution with Chebfun (taking $\varepsilon=0.01$) as follows. We set up the operator and boundary conditions:
N = chebop(-1,1); N.op = @(x,u) 0.01*diff(u,2) + 2*(1-x.^2).*u + u.^2; N.bc = 'dirichlet';
And we set an initial guess:
x = chebfun('x');
N.init = 2*(x.^2-1);
We now solve the problem using the nonlinear backslash operator. By calling \ with two output arguments, we also get the norms of the updates at each iteration returned in a vector. Before solving the problem, we set the cheboppref field 'display' to be 'iter' in order to display information about the solution process. We then plot the solution and the norm of the updates:
cheboppref.setDefaults('display','iter')
[u,info] = solvebvp(N,1);
nrmdu = info.normDelta;
LW = 'linewidth'; MS = 'markersize'; FS = 'fontsize';
subplot(1,2,1), plot(u,LW,1.6), title('Solution',FS,14)
subplot(1,2,2), semilogy(nrmdu,'.-r',LW,1.6,MS,16), title('Convergence',FS,14)
xlim([1 length(nrmdu)]), grid on
Iter. || du || Contraction stepsize len(du) len(u)
--------------------------------------------------------------
01 9.78e-01 6.68e-01 1.0000 64 64
02 2.67e-01 8.70e-02 1.0000 64 64
03 2.85e-02 1.07e-01 1.0000 64 64
04 3.92e-04 1.38e-02 1.0000 64 64
05 7.98e-08 2.03e-04 1.0000 64 64
06 3.57e-15 4.48e-08 1.0000 128 128
--------------------------------------------------------------
Newton's method converged in 6 iterations.
Discretization method used: Collocation.
Final error estimate: 3.57e-15 (differential equation)
0.00e+00 (boundary conditions).
The solution has been obtained to high accuracy:
accuracy = norm(N(u)-1)
accuracy =
8.463254780629571e-14
Running from a different initial guess gives a different solution:
N.init = 2*(x.^2-1).*(1-2./(1+20*x.^2));
cheboppref.setDefaults('display','off')
[u,info] = solvebvp(N,1);
nrmdu = info.normDelta;
subplot(1,2,1), plot(u,LW,1.6), title('Solution',FS,14)
subplot(1,2,2), semilogy(nrmdu,'.-r',LW,1.6,MS,16), title('Convergence',FS,14)
xlim([1 length(nrmdu)]), grid on
accuracy = norm(N(u)-1)
accuracy =
3.126829037542067e-10
2. Alternative boundary conditions
If we want to change the boundary-value problem above to impose different boundary conditions, we can reuse the chebop created and only change the relevant fields. For example, suppose we want to solve the same equation with the new boundary conditions
$$ u(-1) = 1,~~ u'(1) + u(1) = 0. $$
We can execute the following commands:
N.lbc = 1;
N.rbc = @(u) diff(u) + u;
[u,info] = solvebvp(N,1);
nrmdu = info.normDelta;
subplot(1,2,1), plot(u,LW,1.6), title('Solution',FS,14)
subplot(1,2,2), semilogy(nrmdu,'.-r',LW,1.6,MS,16), title('Convergence',FS,14)
xlim([1 length(nrmdu)]), grid on
accuracy = norm(N(u)-1)
accuracy =
3.111051709972451e-10
This example can also be found among the scalar boundary-value problem demos in Chebgui.
References
- C. Bender and S. A. Orzsag, Advanced Mathematical Methods for Scientists and Engineers, McGraw-Hill, 1978 (Section 9.7).