7.1 Introduction
Chebfun has powerful capabilities for solving ordinary differential equations as well as certain partial differential equations. The present chapter is devoted to chebops, the fundamental Chebfun tools for solving ordinary differential (or integral) equations. In particular we focus here on the linear case. We shall see that one can solve a linear two-point boundary value problem to high accuracy by a single backslash command. Nonlinear extensions are described in Section 7.9 and in Chapter 10, and for PDEs in one space dimension, try help pde15s.
Chebfun can also solve certain integral equations, though this topic is not covered much in the Chebfun Guide. See the commands fred and volt and the integro section of the Chebfun Examples collection.
The book Exploring ODEs about ODEs in Chebfun appeared in 2018, and is freely available online as a PDF file [Trefethen, Birkisson & Driscoll 2018]. Appendix B of the book gives 100 short examples of how to solve various ODE problems in Chebfun.
Although one or two examples of initial-value problems for ODEs are presented in this chapter, the emphasis is on boundary-value problems. Beginning with Version 5.1 in December 2014, Chebfun switched to time-stepping methods as the default for initial value problems, a big improvement in speed and robustness in the nonlinear case. See Chapter 10.
7.2 About linear chebops
A chebop represents a differential or integral operator that acts on chebfuns. This chapter focusses on the linear case, though from a user's point of view, linear and nonlinear problems are quite similar. One thing that makes linear operators special is that eigs and expm can be applied to them, as we shall describe in Sections 7.5 and 7.6.
Like chebfuns, chebops start from premise of approximation by piecewise polynomial interpolants; in the context of differential equations, such techniques are called spectral collocation methods. As with chebfuns, the discretizations are chosen automatically to achieve high accuracy. In fact, beginning with version 5, Chebfun actually offers two different methods for solving these problems, which go by the names of rectangular collocation (or Driscoll-Hale) spectral methods and ultraspherical (or Olver-Townsend) spectral methods. See Sections 7.7 and 8.10.
The linear part of the chebop package was conceived at Oxford by Bornemann, Driscoll, and Trefethen [Driscoll, Bornemann & Trefethen 2008], and the original implementation was due to Driscoll, Hale, and Birkisson [Birkisson & Driscoll 2011, Driscoll & Hale 2014]. Much of the functionality of linear chebops is actually implemented in various classes built around the idea of what we call a "linop", but users generally do not deal with linops and related structures directly.
7.3 Chebop syntax
A chebop has a domain, an operator, and sometimes boundary conditions. For example, here is the chebop corresponding to the second-derivative operator on $[-1,1]$:
L = chebop(-1, 1); L.op = @(x,u) diff(u,2);
(For scalar operators like this, one may dispense with the x and just write L.op = @(u) diff(u,2).) This operator can now be applied to chebfuns defined on $[-1,1]$. For example, taking two derivatives of $\sin(3x)$ multiplies its amplitude by 9:
u = chebfun('sin(3*x)');
norm(L(u), inf)
ans = 9.000000000000064
Both the notations L*u and L(u) are allowed, with the same meaning.
min(L*u)
ans = -9.000000000000064
Mathematicians generally prefer to write $Lu$ if $L$ is linear and $L(u)$ if it is nonlinear.
A chebop can also have left and/or right boundary conditions. For a Dirichlet boundary condition it's enough to specify a number:
L.lbc = 0; L.rbc = 1;
More complicated boundary conditions can be specified with anonymous functions, which are then forced to take zero values at the boundary. For example, the following sequence imposes the conditions $u=0$ at the left boundary and $u'=1$ at the right:
L.lbc = @(u) u; L.rbc = @(u) diff(u) - 1;
We can see a summary of L by typing the name without a semicolon:
L
L =
Linear operator:
u |--> diff(u,2)
operating on chebfun objects defined on:
[-1,1]
with
left boundary condition(s):
u = 0
right boundary condition(s):
diff(u)-1 = 0
Boundary conditions are needed for solving differential equations, but they have no effect when a chebop is simply applied to a chebfun. Thus, despite the boundary conditions just specified, L*u gives the same answer as before:
norm(L*u, inf)
ans = 9.000000000000064
When values and derivatives are both to be specified at a single boundary for a scalar ODE, a simpler syntax is also available: instead of writing, say,
L.lbc = @(u) [u-2, diff(u)-3];
one can write
L.lbc = [2; 3];
Here is an example of an integral operator, the operator that maps $u$ defined on $[0,1]$ to its indefinite integral:
L = chebop(0, 1); L.op = @(x,u) cumsum(u);
For example, the indefinite integral of $x$ is $x^2/2$:
x = chebfun('x', [0, 1]);
plot(L*x), grid on
Chebops can be specified in various ways, including all in a single line. For example we could write
L = chebop(@(x,u) diff(u) + diff(u,2), [-1, 1])
L =
Linear operator:
u |--> diff(u)+diff(u,2)
operating on chebfun objects defined on:
[-1,1]
Or we could include boundary conditions:
L = chebop(@(x,u) diff(u) + diff(u,2), [-1, 1], 0, @(u) diff(u))
L =
Linear operator:
u |--> diff(u)+diff(u,2)
operating on chebfun objects defined on:
[-1,1]
with
left boundary condition(s):
u = 0
right boundary condition(s):
diff(u) = 0
For operators applying to more than one variable (needed for solving systems of differential equations), see Section 7.8.
7.4 Solving differential and integral equations
In MATLAB, if A is a square matrix and b is a vector, then the command x=A\b solves the linear system of equations $Ax=b$. Similarly in Chebfun, if L is a differential operator with appropriate boundary conditions and f is a Chebfun, then u=L\f solves the differential equation $L(u)=f$. More generally L might be an integral or integro-differential operator. (Of course, just as you can solve $Ax=b$ only if $A$ is nonsingular, you can solve $L(u)=f$ only if the problem is well-posed.)
For example, suppose we want to solve the differential equation $u''+x^3u = 1$ on the interval $[-3,3]$ with Dirichlet boundary conditions. Here is a Chebfun solution:
L = chebop(-3, 3); L.op = @(x,u) diff(u,2) + x^3*u; L.lbc = 0; L.rbc = 0; u = L\1; plot(u), grid on
We confirm that the computed $u$ satisfies the differential equation to high accuracy:
norm(L(u) - 1)
ans =
1.702191777101623e-11
Let's change the right-hand boundary condition to $u'=0$ and see how this changes the solution:
L.rbc = @(u) diff(u); u = L\1; hold on, plot(u), hold off
An equivalent to backslash is the solvebvp command.
v = solvebvp(L, 1); norm(u - v)
ans =
0
A command like L.bc=100 imposes the corresponding numerical Dirichlet condition at both ends of the domain:
L.bc = 100; plot(L\1), grid on
Boundary conditions can also be specified in a single line, as noted above:
L = chebop( @(x,u) diff(u,2)+10000*u, [-1,1], 0, @(u) diff(u) );
Thus it is possible to set up and solve a differential equation and plot the solution with a single line of Chebfun:
plot( chebop(@(x,u) diff(u,2)+50*(1+sin(x))*u,[-20,20],0,0)\1 ) grid on
When Chebfun solves differential or integral equations, the coefficients may be piecewise smooth rather than globally smooth. For example, here is a 2nd order problem involving a coefficient that jumps from $+1$ (oscillation) for $x<0$ to $-1$ (growth/decay) for $x>0$:
L = chebop(-60, 60); L.op = @(x,u) diff(u,2) - sign(x)*u; L.lbc = 1; L.rbc = 0; u = L\0; plot(u), grid on
Further examples of Chebfun solutions of differential equations with discontinuous coefficients can be found in the Demos menu of chebgui.
Finally, what about periodic boundary conditions? If the boundary condition L.bc='periodic' is specified, Chebfun will discretize the problem by Fourier methods, seeking to find a periodic solution, provided that the right-side function is also periodic. Here is an example:
L = chebop(-pi,pi); L.op = @(x,u) diff(u,2) + diff(u) + 600*(1+sin(x))*u; L.bc = 'periodic'; u = L\1; plot(u), grid on
7.5 Eigenvalue problems: eigs
In MATLAB, eig finds all the eigenvalues of a matrix whereas eigs finds some of them. A differential or integral operator normally has infinitely many eigenvalues, so one could not expect an analog of eig for chebops. eigs, however, has been overloaded. Just like MATLAB eigs, Chebfun eigs finds six eigenvalues by default, together with eigenfunctions if requested. (For details see [Driscoll, Bornemann & Trefethen 2008].) Here is an example involving sine waves.
L = chebop( @(x,u) diff(u,2), [0, pi] ); L.bc = 0; [V, D] = eigs(L); diag(D) clf, plot(V(:,1:4)), ylim([-1 1]), grid on
ans = -35.999999999996703 -25.000000000001538 -16.000000000003045 -8.999999999998460 -4.000000000000856 -0.999999999999610
By default, eigs tries to find the six eigenvalues whose eigenmodes are "most readily converged to", which approximately means the smoothest ones. You can change the number sought and tell eigs where to look for them. Note, however, that you can easily confuse eigs if you ask for something unreasonable, like the largest eigenvalues of a differential operator.
Here we compute 10 eigenvalues of the Mathieu equation and plot the 9th and 10th corresponding eigenfunctions, known as an elliptic cosine and sine. Note the imposition of periodic boundary conditions.
q = 10;
A = chebop(-pi, pi);
A.op = @(x,u) diff(u,2) - 2*q*cos(2*x)*u;
A.bc = 'periodic';
[V, D] = eigs(A, 16, 'LR'); % eigenvalues with largest real part
d = diag(D); [d, ii] = sort(d, 'descend'); V = V(:, ii');
subplot(1,2,1), plot(V(:, 9)), grid on
ylim([-.8 .8]), title('elliptic cosine')
subplot(1,2,2), plot(V(:,10)), grid on
ylim([-.8 .8]), title('elliptic sine')
eigs can also solve generalized eigenproblems, that is, problems of the form $Au = \lambda Bu$. For these one must specify two linear chebops A and B, with the boundary conditions all attached to A. Here is an example of eigenvalues of the Orr-Sommerfeld equation of hydrodynamic linear stability theory at a Reynolds number close to the critical value for eigenvalue instability [Schmid & Henningson 2001]. This is a fourth-order generalized eigenvalue problem, requiring two conditions at each boundary.
Re = 5772;
B = chebop(-1, 1);
B.op = @(x,u) diff(u,2) - u;
A = chebop(-1, 1);
A.op = @(x,u) (diff(u,4) - 2*diff(u, 2) + u)/Re - ...
1i*(2*u + (1 - x^2)*(diff(u, 2) - u));
A.lbc = [0; 0];
A.rbc = [0; 0];
lam = eigs(A, B, 50);
clf, plot(lam, 'r.', 'markersize', 16), grid on, axis equal
spectral_abscissa = max(real(lam))
spectral_abscissa =
-7.806283147737962e-05
For eigenvalue problems of the 1D Schrodinger equation, try help quantumstates.
7.6 Exponential of a linear operator: expm
In MATLAB, expm computes the exponential of a matrix, and this command has been overloaded in Chebfun to compute the exponential of a linear operator. If $L$ is a linear operator and $E(t) = \exp(tL)$, then the partial differential equation $u_t = Lu$ has solution $u(t) = E(t)u(0)$. Thus by taking $L$ to be the 2nd derivative operator, for example, we can use expm to solve the heat equation $u_t = u_{xx}$:
A = chebop(@(x,u) diff(u,2), [-1, 1], 0);
f = chebfun('exp(-1000*(x+0.3)^6)');
clf, plot(f, 'r'), hold on, c = [0.8 0 0];
ylim([-.1 1.1]), grid on
for t = [0.01 0.1 0.5]
u = expm(A, t, f);
plot(u,'color', c), c = 0.5*c;
end
hold off
Here is a more fanciful analogous computation with a complex initial function obtained from the scribble command introduced in Chapter 5.
f = exp(.02i)*scribble('BLUR'); clf
D = chebop(@(x,u) diff(u,2), [-1 1]);
D.bc = 'neumann';
k = 0;
for t = [0 .0001 .001]
k = k + 1; subplot(3,1,k)
if t==0, u = f; else u = expm(D, t, f); end
plot(u, 'linewidth', 3, 'color', [.6 0 1])
xlim(1.05*[-1 1]), axis equal
text(0.01, .46, sprintf('t = %6.4f', t), 'fontsize', 10), axis off
end
7.7 Algorithms: rectangular collocation and ultraspherical
Let us say a word about how Chebfun carries out these computations. Until version 5, Chebfun used classical Chebyshev spectral collocation methods as described in [Trefethen 2000], [Driscoll, Bornemann & Trefethen 2008], and [Driscoll 2010]. With version 5, however, the default changed to a new kind of Chebyshev discretization described in [Aurentz & Trefethen 2017], [Driscoll & Hale 2014] and [Xu & Hale 2014]. These rectangular collocation or Driscoll-Hale spectral discretizations start from the idea that a differential operator is discretized as a rectangular matrix that maps from one grid to another with fewer points. The matrix is then made square again by the incorporation of boundary conditions. When a differential equation is solved in Chebfun, the problem is solved on a sequence of grids of size $33, 65, \dots$ until convergence is achieved in the usual Chebfun sense defined by decay of Chebyshev expansion coefficients.
If you want to learn the details of rectangular discretizations, you can find a sequence of 12 explicit Chebfun examples presented in [Aurentz & Trefethen 2017]. As described there, you can get your hands on Chebfun's discretization matrices with the command matrix. For example, here is the $6\times 6$ discretization matrix for the second derivative operator on $[-1,1]$ with zero boundary conditions:
L = chebop(@(u) diff(u,2)); L.bc = 0; format short matrix(L,4) format long
ans =
1.0000 0 0 0 0 0
0 0 0 0 0 1.0000
32.5709 -51.8814 27.9770 -14.2643 10.1922 -4.5944
-0.5502 5.5850 -10.1042 7.0354 -3.3398 1.3737
1.3737 -3.3398 7.0354 -10.1042 5.5850 -0.5502
-4.5944 10.1922 -14.2643 27.9770 -51.8814 32.5709
The first two rows correspond to the boundary conditions, and the remaining $4\times 6$ block is the rectangular discretization matrix that takes input from the 6-point Chebyshev grid, interpolates it by a degree-5 polynomial, differentiates the interpolant twice, and samples the result on the 4-point Chebyshev grid.
One matter you might not guess was challenging is the determination of whether or not an operator is linear! This is important since if an operator is linear, special actions are possible possible such as application of eigs and expm and solution of differential equations in a single step without iteration. Chebfun includes special devices to determine whether a chebop is linear so that these effects can be realized [Birkisson 2014].
As mentioned, the discretization length of a Chebfun solution is chosen automatically according to the intrinsic resolution requirements. However, the matrices that arise in Chebyshev spectral methods are notoriously ill-conditioned. Thus the final accuracy in solving differential equations in Chebfun's default mode is rarely close to machine precision. Typically one loses two or three digits for second-order differential equations and five or six for fourth-order problems.
This problem of ill-conditioning was one of the motivations for the development of the other discretization method used by Chebfun, known as ultraspherical or Olver-Townsend spectral discretizations [Olver & Townsend 2013]. Here, rather than a single Chebyshev basis, several different bases of ultraspherical polynomials are used, depending on the order of the differential operator. This leads to better conditioned matrices that are also sparser, especially for linear problems with constant or smooth coefficients. By default, Chebfun uses rectangular collocation discretizations, but most problems can equally be solved in ultraspherical mode, and the results will sometimes be more accurate. For example, here we solve a problem whose exact solution is $\cos(x)$ in the rectangular fashion and check the error at $x=5$:
tic u = chebop(@(x, u) diff(u, 2) + u, [-10,10], cos(10), cos(10))\0; toc error = u(5) - cos(5)
Elapsed time is 0.068711 seconds.
error =
-2.570166302007237e-14
We can switch to ultraspherical mode and run the same experiment again like this:
cheboppref.setDefaults('discretization', @ultraS)
tic
u = chebop(@(x,u) diff(u,2) + u, [-10,10], cos(10), cos(10))\0;
toc
error = u(5) - cos(5)
cheboppref.setDefaults('factory') % reset to standard mode
Elapsed time is 0.281529 seconds.
error =
-1.387778780781446e-15
7.8 Block operators and systems of equations
Some problems involve several variables coupled together. In Chebfun, these are treated with the use of quasimatrices, that is, chebfuns with several columns, as described in Chapter 6.
For example, suppose we want to solve the coupled system $u'=v$, $v'=-u$ with initial data $u=1$ and $v=0$ on the interval $[0,10\pi]$. (This comes from writing the equation $u''=-u$ in first-order form, with $v=u'$.) We can solve the problem like this:
L = chebop(0, 10*pi); L.op = @(x, u, v) [diff(u) - v; diff(v) + u]; L.lbc = @(u, v) [u-1; v]; rhs = [0; 0]; U = L\rhs;
The solution U is an $\infty\times 2$ Chebfun quasimatrix with columns u=U(:,1) and v=U(:,2). Here is a plot:
clf, plot(U), grid on
The overloaded spy command helps clarify the structure of the operator we just made use of:
spy(L)
This image shows that $L$ maps a pair of functions $[u;v]$ to a pair of functions $[w;y]$, where the dependences of $w$ on $u$ and $y$ on $v$ are global (because of the derivative) whereas the dependences of $w$ on $v$ and $y$ on $u$ are local (diagonal).
To illustrate the solution of an eigenvalue problem involving a block operator, we can take much the same idea. The eigenvalue problem $u''=c^2u$ with $u=0$ at the boundaries can be written in first order form as $u'=cv$, $v'=cu$.
L = chebop(0, 10*pi); L.op = @(x, u, v) [diff(v); diff(u)]; L.lbc = @(u,v) u; L.rbc = @(u,v) u;
The operator in this eigenvalue problem has a simpler structure than before:
clf, spy(L)
Here are the first 7 eigenvalues:
[eigenfunctions,D] = eigs(L, 7); eigenvalues = diag(D)
eigenvalues = -0.000000000000001 + 0.000000000000000i -0.000000000000000 - 0.100000000000000i -0.000000000000000 + 0.100000000000000i 0.000000000000001 - 0.200000000000001i 0.000000000000001 + 0.200000000000001i -0.000000000000000 - 0.300000000000002i -0.000000000000000 + 0.300000000000002i
The eigenfunctions result has the first seven eigenfunctions for each of the two variables, u and v. We could extract this chebmatrix result to a chebfun like this:
U = chebfun( eigenfunctions(1,:) ); V = chebfun( eigenfunctions(2,:) ); size(V)
ans = Inf 7
Both U and V are complex, but only because of roundoff:
normRealU = norm(real(U)) normImagV = norm(imag(V))
normRealU =
7.203750542739416e-12
normImagV =
7.210253604889195e-12
This fact makes it easy to plot them.
subplot(2,1,1)
plot(imag(U)), ylabel('imag(U)')
subplot(2,1,2)
plot(real(V)), ylabel('real(V)')
7.9 Nonlinear equations by Newton iteration
As mentioned at the beginning of this chapter, nonlinear differential equations are discussed in Chapter 10. As an indication of some of the possibilities, however, we now illustrate how a sequence of linear problems may be useful in solving nonlinear problems. For example, the nonlinear BVP $$ 0.001u'' - u^3 = 0,\qquad u(-1)=1,~~ u(1)=-1 $$ could be solved by Newton iteration as follows.
L = chebop(-1, 1);
L.op = @(x,u) 0.001*diff(u, 2);
J = chebop(-1, 1);
x = chebfun('x');
u = -x; nrmdu = Inf;
while nrmdu > 1e-10
r = L*u - u.^3;
J.op = @(du) .001*diff(du, 2) - 3*u^2*du;
J.bc = 0;
du = -(J\r);
u = u + du; nrmdu = norm(du)
end
clf, plot(u), grid on
nrmdu =
0.260668532007020
nrmdu =
0.164126069559937
nrmdu =
0.098900892365439
nrmdu =
0.053787171683933
nrmdu =
0.021518152858429
nrmdu =
0.003586696693249
nrmdu =
8.951602489054714e-05
nrmdu =
5.357404836265776e-08
nrmdu =
1.709262919444643e-14
Note the beautifully fast convergence, as one expects with Newton's method. The chebop J defined in the while loop is a Jacobian operator (=Frechet derivative), which we have constructed explicitly by differentiating the nonlinear operator defining the ODE. In Section 10.4 we shall see that this whole Newton iteration can be automated by use of Chebfun's "nonlinear backslash" capability, which utilizes automatic differentiation to construct the Frechet derivative automatically. In fact, all you need to type is
N = chebop(-1, 1); N.op = @(x,u) 0.001*diff(u, 2) - u^3; N.lbc = 1; N.rbc = -1; v = N\0;
The result is the same as before to many digits of accuracy:
norm(u - v)
ans =
1.201029224139150e-12
7.10 BVP systems with unknown parameters
Sometimes ODEs or systems of ODEs contain unknown parameter values that must be computed as part of the solution. An example of this is MATLAB's built-in mat4bvp example. These parameters can always be included in a system as unknowns with zero derivatives, but this can be computationally inefficient. Chebfun allows the option of explicit treatment of the parameters. Often the dependence of the solution on these parameters is nonlinear (as in the case below), and this discussion might better have been left to Chapter 10, but since, from the user perspective, there is little difference in this case, we include it here.
Below is an example of such a parameterised problem, which represents a linear pendulum with a forcing sine-wave term of an unknown frequency $T$. The task is to compute the solution for which $$ u(-\pi) = u(\pi) = u'(\pi) = 1. $$
N = chebop(@(x, u , T) diff(u,2) - u - sin(T*x/pi), [-pi pi]); N.lbc = @(u,T) u - 1; N.rbc = @(u,T) [u - 1; diff(u) - 1]; uT = N\0;
Here, the output uT is a chebmatrix -- an object that is amongst other features able to vertically concatenate chebfuns and scalar. The first entry corresponds to the variable $u$ while the second is the scalar $T$. We could access the elements of uT via curly-brackets syntax,
u = uT{1}; T = uT{2};
but one can access the same results more simply like this:
[u,T] = N\0; T plot(u), grid on
T = 0.005438812795290
As the system is nonlinear in $T$, we can expect that there will be more than one solution. Indeed, if we choose a different initial guess for $T$, we can converge to one of these.
N.init = [chebfun(1, [-pi pi]); 4]; [u,T] = N\0; T = T(1) plot(u), grid on
T = 4.044049959218974
7.11 References
[Aurentz & Trefethen 2017] J. L. Aurentz and L. N. Trefethen, Block opearators and spectral discretizations, SIAM Review 59 (2017), 423--446.
[Birkisson 2014] A. Birkisson, Numerical Solution of Nonlinear Boundary Value Problems for Ordinary Differential Equations in the Continuous Framework, D. Phil. thesis, University of Oxford, 2014.
[Birkisson & Driscoll 2011] A. Birkisson and T. A. Driscoll, "Automatic Frechet differentiation for the numerical solution of boundary-value problems", ACM Transactions on Mathematical Software, 38 (2012), 1-26.
[Driscoll 2010] T. A. Driscoll, "Automatic spectral collocation for integral, integro-differential, and integrally reformulated differential equations", Journal of Computational Physics, 229 (2010), 5980-5998.
[Driscoll, Bornemann & Trefethen 2008] T. A. Driscoll, F. Bornemann, and L. N. Trefethen, "The chebop system for automatic solution of differential equations", BIT Numerical Mathematics, 46 (2008), 701-723.
[Driscoll & Hale 2014] T. A. Driscoll and N. Hale, "Rectangular spectral collocation", IMA Journal of Numerical Analysis 36 (2016), 108--132.
[Fornberg 1996] B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge University Press, 1996.
[Olver & Townsend 2013] S. Olver and A. Townsend, "A fast and well-conditioned spectral method," SIAM Review, 55 (2013), 462-489.
[Schmid & Henningson 2001] P. J. Schmid and D. S. Henningson, Stability and Transition in Shear Flows, Springer, 2001.
[Trefethen 2000] L. N. Trefethen, Spectral Methods in MATLAB, SIAM, 2000.
[Trefethen, Birkisson & Driscoll 2018] L. N. Trefethen, A. Birkisson, and T. A. Driscoll, Exploring ODEs, SIAM, 2018. Freely available at http://people.maths.ox.ac.uk/trefethen/ExplODE/.
[Xu & Hale 2015] K. Xu and N. Hale, "Explicit construction of rectangular differentiation matrices", IMA Journal of Numerical Analysis 36 (2016), 618-632.