2.2.1. Approximation of the solution by collocation

For engineers, a practical question concerns the determination of approximate solutions. For instance, the solution may be determined at some particular values of t – for instance, t₁,…,t_ns – leading to particular values of y− for instance, x₁,…,x_ns – and an interpolation may be performed in order to generate values for arbitrary t. This approach is usually referred to as collocation. In practice, it involves a family of interpolation functions

such as, for instance, φ_i(t) = tⁱ⁻¹ (polynomial interpolation). Then, an approximation:

is determined by solving the linear system of algebraic equations:

The unknowns to be determined are the coefficients u_p ∈ ⁿ, for p = 1, … k. They can be arranged in a matrix such that u_pj = (u_p)_j However, the data x_r ∈ ⁿ may be arranged in a matrix such that x_ri = (x_r)_i, so that these equations read as:

where, for 1 ≤ i,j ≤ n, 1≤ p ≤ k, 1≤ r ≤ ns,

These equations may be transformed in a standard linear system involving two indexes by using the transformation.

By setting I = index(r, i, ns), J = index(p, j, k), , , , we have . This linear system has dimension n.ns × k.n. For ns > k, it is overdetermined and least squares solution may be used – under Matlab^®, this is made automatically. Notice that it may also be solved component by component: let us introduce:

Then, for 1 ≤ i ≤ n,

Thus, the solution may be sequentially determined by solving n linear systems of size ns × k.

2.2.2. Variational approximation of the solution

Collocation may be sensitive to errors in the numerical determination of x₁, … x_ns. The variational formulations presented furnish a more robust method for the determination of the solution.

By denoting Φ = (φ₁, … ,φ_k)^t, we have x ≈ U^tΦ. By taking y_i(t) = φ_j(t) in equation [2.6], we obtain, for 1 ≤ i ≤ m; 1 ≤ p ≤ k.

[2.7]

so that;

[2.8]

This system contains k × m equations for k × n unknowns. When m = n, the numbers of equations and unknowns coincide For m < n, the equations are under-determined (less equations than unknowns); for m > n, the equations are over-determined (more equations than unknowns). As in the previous situation, this classification does not imply the existence or nonexistence of solutions. Equation [2.7] must be solved by an appropriated method, such as, for instance, Newton–Raphson, quasi-newton or fixed-point iterations. As an alternative, we may look for:

[2.9]

We observe that;

[2.10]

with

[2.11]

These equations may be invoked for the evaluation of the gradients when using optimization methods.

2.2.3. Linear equations

The particular situation where

leads to a linear system. In this case,

so that;

and

where

These equations may be transformed in a standard linear system involving two indexes by using the transformation previously introduced: I = index(p, i, k), J = index(q, j, k), , , . Then, .

2.2.4. Connection to orthogonal projections

The variational approach may be interpreted in terms of orthogonal projection: let V be the set of all the functions defined on (a, b) and taking their values on .

V is a linear space (i.e. a vector space) and generates a finite dimensional subspace S ⊂ V:

Since dim(S) = k, S may be assimilated to ^k: the element s = s^tΦ∈ S is entirely defined by the vector s = (s₁,…,s_k)^t ∈ ^k. Let us consider a particular scalar product for vectors u, v ∈ ^k:

We have:

where = u^tΦ, v = v^tΦ ∈ S ⊂ V. Let us extend this definition as being a scalar product on the whole space V:

Recall that u ⊥ v ⇔ (u, v) = 0, i.e. u and v; are orthogonal if and only if their scalar product is null. Let P_S:V → S denote the orthogonal projection onto S. For u ∈ V, P_s(u) is characterized by (see section 3.7)

or, equivalently

Equation [2.7] shows that, for 1 ≤ i ≤ m,

So, for 1≤ i ≤ m,

Thus,

So, the orthogonal projection of F_i(U^tΦ(t), t) onto S is null for 1 ≤ i ≤ m:

Let us consider the Cartesian product S^m = S × … × S (m times) We have also

Thus, the orthogonal projection of F(U^tΦ(t), t) onto S^m is null:

In the particular situation where

we have, on the one hand, m = n and, on the other hand,

Thus, and, on the one hand,

is the orthogonal projection of f_i onto S

while, on the other hand (recall that m = n),

is the orthogonal projection of f onto S^m.

2.2.5. Numerical determination of the orthogonal projections

Let us consider an element u ∈ Vⁿ and its orthogonal projection onto Sⁿ We have:

Thus, ,where;

So,

Analogously to collocation, these equations form a linear system for the coefficients . By setting, for 1 ≤ p, q ≤ k, 1 ≤ i, j ≤ n,

we have:

These equations may be transformed in a standard linear system involving two indexes by setting I = index(q, i, k),

J = index(p, j, k), , , : then . This linear system has dimension kn × kn. It may also be solved component by component: let us introduce:

Then, for 1 ≤ i ≤ n,

Thus, the solution may be sequentially determined by solving n linear systems of size × k.

2.2.6. Matlab^® classes for a numerical solution

In order to determine numerical solutions, we need to define the family For instance, in one-dimensional (1D) situations, we may use a polynomial family defined by the following class:

In this code, the polynomials are , with xmin=a and xmin= b. Property zero defines numbers which are treated as zero. degree defines the degree of the polynomial to be used. Property spmatrix is a square matrix containing the integrals of the products of the elements of the basis, while property intmatrix is a vector containing the integrals of the elements themselves.

Let us assume that the data is given by a structure approxdata having as properties: points, values, subprogram, which contain, respectively, the points x, the corresponding values f(x) and a subprogram evaluating f(x). In this case, the coefficients of an approximation and, namely, orthogonal projections may be evaluated by using the following class:

This class contains a method approx_coefficients which determines the coefficients of the approximation according to the method chosen: ‘collocation’ uses the collocation approach, while the others use the variational approach, but differ in the manner where the integrals are evaluated: ‘variational_mean’ uses the data in order to evaluate a mean value approximation, ‘variational_int’ uses the quadrature by Matlab®, ‘variational_sp’ uses the matrix of integrals spmatrix furnished in the definition of the polynomial basis. The class also contains methods integrate and derive for the evaluation of an integral on the interval (xmin, xmax) and the derivative at a vector of points, respectively.

Let us illustrate the use of these classes. For instance, let us consider the noisy data generated as follows:

f = @(t) exp(t);
t = (-1:0.2:1);
p.x = t';
p.dim = 1;
x = spam.partition(f,p);
noise = 0.1;
td = (-1:0.05:1);
xd = spam.points('vector',f,td);
vd = (1+ noise*(2*rand(size(xd))-1)).*xd;
fnoise = @(t) interp1(td,vd,t,'linear');
data.points = p;
data.values = (1+ noise*(2*rand(size(x))-1)).*x;
data.subprogram = @(t) fnoise(t);
data.dim = 1;

These commands generate noisy data representing the exponential function on the interval (−1,1). The level of noise is 10%. Notice that the subprogram evaluating the function is not the exact exponential, but a noisy version generated by linear interpolation of noisy data. The commands

xmin = -1;
xmax = 1;
degree = 4;
zero = 1e-3;
phi = polynomial_basis(zero,xmin,xmax,degree);
b = basis(phi);
u1 = b.approx_coeffs('collocation',data);

generate the coefficients of the approximation by collocation by using a polynomial basis involving a maximal degree 4. The command

generates the coefficients corresponding to the variational approximation and evaluate the integrals by using the means. The command

generates the variational approximation and evaluates the integrals by using subprogram integration.

The results are shown in Figure 2.2 (recall that data is given on (−1,1) and points outside this interval are extrapolated ones):

The derivative is shown below.

We observe that, on the one hand, the results by collocation and variational approach with integrals evaluated by the mean are close and, on the other hand, the results furnished by the variational approach with more sophisticated integration are close to exact ones. Concerning the evaluation of the integral of f on the interval (−1,1), the commands

v1 = b.integrate(u1);
v2 = b.integrate(u2);
v3 = b.integrate(u3);
v4 = b.integrate(u4);
ff = @(t) b.projection(u4,t);
v5 = integral(f,-1,1);

furnish v1 = v2 = 2.3535, v3 = v4 = v5 = 2.3442, while the exact value is e − e⁻¹≈ 2.3504.

Results for the vector (t) = (sin(πt) exp(t)^t are given in Figures 2.4 and 2.5. We use a degree 5. In this situation, we obtain v1 = v2 = (−0.0139 2.3507), v3 = 0.0027 2.3658).

When considering nonlinear parametric equations involving a single parameter, we may use the following class

In this class, the solution is obtained either by minimizing the quadratic mean square of the elements of (U) or by solving the equations (U) = 0. The optimization step is performed by using the intrinsic Matlab® function fminsearch., while the determination of a zero is performed by the intrinsic function fzero. Parameters such as number of iterations, function evaluations, precision are defined by using optimset. Choice ‘mean’ evaluates the integrals by using the mean value approximation with ns equal subintervals, while ‘integral’ uses the intrinsic quadrature of Matlab®. The equations are given by structured data in equations, which has as fields subprogram, neq, dim, which correspond, respectively, to the subprogram which furnishes the value of F(x, t), the number of equations and the number of unknowns (dimension of the unknown).

EXAMPLE 2.3.–Let us consider the simple situation where,

We use the polynomial basis with a degree 6. The code is:

xmin = 0;
xmax = pi;
ns = 500;
degree = 6;
zero = 1e-3;
phi = polynomial_basis(zero,xmin,xmax,degree);
b = basis(phi);
f = @(x,t) x - sin(t);
eqs.subprogram = f;
eqs.dim = 1;
eqs.neq = 1;

The obvious solution is x = sin(t). The results furnished by using:

pa = parametric_algebraic(eqs,b,ns);
nitmax = 10000;
U1 = pa.solve_zero('integral',U0,nitmax);
U2 = pa.solve_min('mean',U0,nitmax)

are shown in Figure 2.6 below (left: methods integral and solve_zero, right: methods mean and solve_min).

EXAMPLE 2.4.– Let us consider the simple situation where,

We use the polynomial basis with a degree 5. The code is:

xmin = 0;
xmax = pi;
ns = 500;
degree = 5;
zero = 1e-3;
phi = polynomial_basis(zero,xmin,xmax,degree);
b = basis(phi);
f = @(x,t) [x(1) + x(2) - sin(t) ; x(1) - x(2) -
cos(t)];
eqs.subprogram = f;
eqs.dim = 2;
eqs.neq = 2;
pa = parametric_algebraic(eqs,b,ns);
U1 = pa.solve_zero(
'mean',U0,nitmax);
U2 = pa.solve_min(
'integral',U0,nitmax);

The obvious solution is x₁ = (sin(t) + cos (t))/2, x₁ = (sin(t) − cos(t))/2. The results are shown in Figures 2.7 (using method solve_zero, integrals evaluated by the mean) and 2.8 (using method solve_min, integrals evaluated by quadrature) below.

EXAMPLE 2.5.– Let us consider the simple situation where

We use the polynomial basis with a degree 6. The results furnished by using method solve_zero are shown in Figure 2.9.

EXAMPLE 2.6.– Let us consider the simple situation where

We use the polynomial basis with a degree 5. The results are shown in Figure 2.10 (using method solve_zero).

Other families may be used, such as, for instance, the trigonometric family: φ₁(t) = 1, φ_2i(t) = cos(it) φ_2i+1(t) = sin(it) or the family associated with finite elements P1. These families are defined by the classes below.

Let us illustrate the use of these classes: consider again f(x, t) = x−sin(t), what corresponds to the orthogonal projection of x = sin (t) onto the linear subspaces generated by the basis. The results are shown in the figures below. P1 Finite Elements have a piecewise constant derivative.