All the concepts and algorithms developed for approximation of 1D functions f(x) can readily be extended to 2D functions f(x,y) and 3D functions f(x,y,z). Basically, the extensions consist of defining basis functions ψi(x,y) or ψi(x,y,z) over some domain Ω, and for the least squares and Galerkin methods, the integration is done over Ω.
As in 1D, the least squares and projection/Galerkin methods lead to linear systems ∑j∈IsAi,jcj=bi,i∈Is,Ai,j=(ψi,ψj),bi=(f,ψi), where the inner product of two functions f(x,y) and g(x,y) is defined completely analogously to the 1D case (24): (f,g)=∫Ωf(x,y)g(x,y)dxdy
One straightforward way to construct a basis in 2D is to combine 1D basis functions. Say we have the 1D vector space Vx=span{ˆψ0(x),…,ˆψNx(x)}. A similar space for a function's variation in y can be defined, Vy=span{ˆψ0(y),…,ˆψNy(y)}. We can then form 2D basis functions as tensor products of 1D basis functions.
Tensors are typically represented by arrays in computer code. In the above example, a and b are represented by one-dimensional arrays of length M and N, respectively, while p=a⊗b must be represented by a two-dimensional array of size M×N.
Tensor products can be used in a variety of context.
Given the vector spaces Vx and Vy as defined in (115) and (116), the tensor product space V=Vx⊗Vy has a basis formed as the tensor product of the basis for Vx and Vy. That is, if {φi(x)}i∈Ix and {φi(y)}i∈Iy are basis for Vx and Vy, respectively, the elements in the basis for V arise from the tensor product: {φi(x)φj(y)}i∈Ix,j∈Iy. The index sets are Ix={0,…,Nx} and Iy={0,…,Ny}.
The notation for a basis function in 2D can employ a double index as in ψp,q(x,y)=ˆψp(x)ˆψq(y),p∈Ix,q∈Iy. The expansion for u is then written as a double sum u=∑p∈Ix∑q∈Iycp,qψp,q(x,y). Alternatively, we may employ a single index, ψi(x,y)=ˆψp(x)ˆψq(y), and use the standard form for u, u=∑j∈Iscjψj(x,y). The single index is related to the double index through i=p(Ny+1)+q or i=q(Nx+1)+p.
Suppose we choose ˆψp(x)=xp, and try an approximation with Nx=Ny=1: ψ0,0=1,ψ1,0=x,ψ0,1=y,ψ1,1=xy. Using a mapping to one index like i=q(Nx+1)+p, we get ψ0=1,ψ1=x,ψ2=y,ψ3=xy.
With the specific choice f(x,y)=(1+x2)(1+2y2) on Ω=[0,Lx]×[0,Ly], we can perform actual calculations: A0,0=(ψ0,ψ0)=∫Ly0∫Lx0ψ0(x,y)2dxdy=∫Ly0∫Lx0dxdy=LxLy,A1,0=(ψ1,ψ0)=∫Ly0∫Lx0xdxdy=12L2xLy,A0,1=(ψ0,ψ1)=∫Ly0∫Lx0ydxdy=12L2yLx,A0,1=(ψ0,ψ1)=∫Ly0∫Lx0xydxdy=∫Ly0ydy∫Lx0xdx=14L2yL2x. The right-hand side vector has the entries b0=(ψ0,f)=∫Ly0∫Lx01⋅(1+x2)(1+2y2)dxdy=∫Ly0(1+2y2)dy∫Lx0(1+x2)dx=(Ly+23L3y)(Lx+13L3x)b1=(ψ1,f)=∫Ly0∫Lx0x(1+x2)(1+2y2)dxdy=∫Ly0(1+2y2)dy∫Lx0x(1+x2)dx=(Ly+23L3y)(12L2x+14L4x)b2=(ψ2,f)=∫Ly0∫Lx0y(1+x2)(1+2y2)dxdy=∫Ly0y(1+2y2)dy∫Lx0(1+x2)dx=(12Ly+12L4y)(Lx+13L3x)b3=(ψ2,f)=∫Ly0∫Lx0xy(1+x2)(1+2y2)dxdy=∫Ly0y(1+2y2)dy∫Lx0x(1+x2)dx=(12L2y+12L4y)(12L2x+14L4x).
There is a general pattern in these calculations that we can explore. An arbitrary matrix entry has the formula Ai,j=(ψi,ψj)=∫Ly0∫Lx0ψiψjdxdy=∫Ly0∫Lx0ψp,qψr,sdxdy=∫Ly0∫Lx0ˆψp(x)ˆψq(y)ˆψr(x)ˆψs(y)dxdy=∫Ly0ˆψq(y)ˆψs(y)dy∫Lx0ˆψp(x)ˆψr(x)dx=ˆA(x)p,rˆA(y)q,s, where ˆA(x)p,r=∫Lx0ˆψp(x)ˆψr(x)dx,ˆA(y)q,s=∫Ly0ˆψq(y)ˆψs(y)dy, are matrix entries for one-dimensional approximations. Moreover, i=qNy+q and j=sNy+r.
With ˆψp(x)=xp we have ˆA(x)p,r=1p+r+1Lp+r+1x,ˆA(y)q,s=1q+s+1Lq+s+1y, and Ai,j=ˆA(x)p,rˆA(y)q,s=1p+r+1Lp+r+1x1q+s+1Lq+s+1y, for p,r∈Ix and q,s∈Iy.
Corresponding reasoning for the right-hand side leads to bi=(ψi,f)=∫Ly0∫Lx0ψifdxdx=∫Ly0∫Lx0ˆψp(x)ˆψq(y)fdxdx=∫Ly0ˆψq(y)(1+2y2)dy∫Ly0ˆψp(x)xp(1+x2)dx=∫Ly0yq(1+2y2)dy∫Ly0xp(1+x2)dx=(1q+1Lq+1y+2q+3Lq+3y)(1p+1Lp+1x+2q+3Lp+3x)
Choosing Lx=Ly=2, we have A=[444441634163441631634163163649],b=[308914034460],c=[−1943−238]. Figure 35 illustrates the result.
Figure 35: Approximation of a 2D quadratic function (left) by a 2D bilinear function (right) using the Galerkin or least squares method.
The least_squares
function from
the section Orthogonal basis functions and/or the
file approx1D.py
can with very small modifications solve 2D approximation problems.
First, let Omega
now be a list of the intervals in x and y direction.
For example, Ω=[0,Lx]×[0,Ly] can be represented
by Omega = [[0, L_x], [0, L_y]]
.
Second, the symbolic integration must be extended to 2D:
import sympy as sym
integrand = psi[i]*psi[j]
I = sym.integrate(integrand,
(x, Omega[0][0], Omega[0][1]),
(y, Omega[1][0], Omega[1][1]))
provided integrand
is an expression involving the sympy
symbols x
and y
.
The 2D version of numerical integration becomes
if isinstance(I, sym.Integral):
integrand = sym.lambdify([x,y], integrand)
I = sym.mpmath.quad(integrand,
[Omega[0][0], Omega[0][1]],
[Omega[1][0], Omega[1][1]])
The right-hand side integrals are modified in a similar way.
(We should add that sympy.mpmath.quad
is sufficiently fast
even in 2D, but scipy.integrate.nquad
is much faster.)
Third, we must construct a list of 2D basis functions. Here are two examples based on tensor products of 1D "Taylor-style" polynomials xi and 1D sine functions sin((i+1)πx):
def taylor(x, y, Nx, Ny):
return [x**i*y**j for i in range(Nx+1) for j in range(Ny+1)]
def sines(x, y, Nx, Ny):
return [sym.sin(sym.pi*(i+1)*x)*sym.sin(sym.pi*(j+1)*y)
for i in range(Nx+1) for j in range(Ny+1)]
The complete code appears in approx2D.py.
The previous hand calculation where a quadratic f was approximated by a bilinear function can be computed symbolically by
>>> from approx2D import *
>>> f = (1+x**2)*(1+2*y**2)
>>> psi = taylor(x, y, 1, 1)
>>> Omega = [[0, 2], [0, 2]]
>>> u, c = least_squares(f, psi, Omega)
>>> print u
8*x*y - 2*x/3 + 4*y/3 - 1/9
>>> print sym.expand(f)
2*x**2*y**2 + x**2 + 2*y**2 + 1
We may continue with adding higher powers to the basis:
>>> psi = taylor(x, y, 2, 2)
>>> u, c = least_squares(f, psi, Omega)
>>> print u
2*x**2*y**2 + x**2 + 2*y**2 + 1
>>> print u-f
0
For Nx≥2 and Ny≥2 we recover the exact function f, as expected, since in that case f∈V (see the section Perfect approximation).
Extension to 3D is in principle straightforward once the 2D extension is understood. The only major difference is that we need the repeated outer tensor product, V=Vx⊗Vy⊗Vz. In general, given vectors (first-order tensors) a(q)=(a(q)0,…,a(q)Nq, q=0,…,m, the tensor product p=a(0)⊗⋯⊗am has elements pi0,i1,…,im=a(0)i1a(1)i1⋯a(m)im. The basis functions in 3D are then ψp,q,r(x,y,z)=ˆψp(x)ˆψq(y)ˆψr(z), with p∈Ix, q∈Iy, r∈Iz. The expansion of u becomes u(x,y,z)=∑p∈Ix∑q∈Iy∑r∈Izcp,q,rψp,q,r(x,y,z). A single index can be introduced also here, e.g., i=NxNyr+qNx+p, u=∑iciψi(x,y,z).