Lecture 4: The weighted residual and weak formulations

We have seen that the method of weighted residuals consists of taking our differential equation:
				A(u) = f(x)  in Ω 
where Ω is the domain where the differential equation holds, and requiring the residual that results from applying the differential equation to our approximate solution (Subscript[U, N]): 

	(Subscript[U, N])= (∑_(j=1)^N)(Subscript[c, j])(Subscript[ϕ, j])(x)+(Subscript[ϕ, 0])(x)
	R((Subscript[U, N])) =A((Subscript[U, N])) -f= A((∑_(j=1)^N)(Subscript[c, j])(Subscript[ϕ, j])(x)+(Subscript[ϕ, 0])(x)) - f
to vanish in an weighted integral sense:
				(∫_Ω)w R((Subscript[U, N]))dΩ=0

If A(u) is a linear operator, i.e., it depends linearly on "u", then we have:

			R((Subscript[U, N]))=A((Subscript[U, N])) = (Underoverscript[Σ, (j=1), N])(Subscript[c, j])A((Subscript[ϕ, j])) - f
(Subscript[U, N]) has "N" unknown coefficients (Subscript[c, i]), i=1,N. In order to determine them, we need to compute "N" integrals with "N" different "w" functions. 
		(∫_Ω)(Subscript[w, i]) R((Subscript[U, N]))dΩ= (∫_Ω)(Subscript[w, i])((Underoverscript[Σ, (j=1), N])(Subscript[c, j])A((Subscript[ϕ, j])) - f)dΩ = 
(Underoverscript[Σ, (j=1), N])(Subscript[c, j])[(∫_Ω)(Subscript[w, i]) A((Subscript[ϕ, j]))dΩ]-(∫_Ω)(Subscript[w, i]) fdΩ = 0, i=1,N
After the integrals are evaluated, a system of algebraic equations on the coefficients are obtained, from which the approximate solution (Subscript[U, N]) can be obtained. The expression inbetween the square brackets is the so-called "stiffness matrix"

			(Subscript[K, ij])= (∫_Ω)(Subscript[w, i]) A((Subscript[ϕ, j]))dΩ, i,j=1,N
In general, the stiffness matrix is not symmetric. The second integral gives the "force" array whose components are:
			(Subscript[r, i])=(∫_Ω)(Subscript[w, i]) fdΩ  

Then, the system above can be written in matrix form as:

			[K] {U} = {R}
where
(({U})^T)= {(Subscript[c, 1]) (Subscript[c, 2]) ..., (Subscript[c, N])}
(({R})^T)= {(Subscript[r, 1])(Subscript[r, 2])..., (Subscript[r, N])}

Note that in this case (Subscript[U, N])has to:
	- be smooth enough such that all the derivatives appearing in A() can be computed
	- satisfy the boundary conditions of the problem
There are much less stringent requirements on "w": No derivatives, no boundary condition requirements.
The weighted-residual statement "replaces" the differential equation. The boundary conditions need to be applied separately by properly choosing the approximation function (Subscript[U, N]).

Example:
		A() = ((-d)/dx)(a(x) (d/dx)());  Ω ≡{x ∈ ℝ, 0<x<1}
		((-d)/dx)(a(x) (du/dx))=f(x); 0<x<1; u(0) = (Subscript[u, 0]); a(x) (du/dx)(Subscript[|, x=L])= (Subscript[Q, 0])
		
		weighted residual formulation:
		
		(∫_0^1)w (((-d)/dx)(a(x) ((Subscript[dU, N])/dx))-f(x))dx = 0

How do we choose the weight functions?
There are various ways, each giving rise to a different method:

Galerkin method
the weight functions are chosen to be the same as the approximation functions:

		w(x) = (Subscript[ϕ, i])(x), i=1,N, replace in weighted residual integral:
		
		(Underoverscript[Σ, (j=1), N])(Subscript[c, j])[(∫_Ω)(Subscript[ϕ, i]) A((Subscript[ϕ, j]))dΩ]-(∫_Ω)(Subscript[ϕ, i]) fdΩ = 0, i=1,N

Note: the stiffness matrix is not symmetric here either.
Example:
	((-(d^2)u)/(dx^2))- u + (x^2)=0, u(0) = 0, u'(1)=1
	For a weighted residual method, the interpolation functions should satisfy:
	(Subscript[ϕ, 0])(0) = 0, (Subscript[ϕ, 0])'(1) = 1 require (Subscript[ϕ, 0]) to satisfy all the actual boundary conditions
	(Subscript[ϕ, j])(0) = 0, (Subscript[ϕ, j])'(1) = 0 satisfy homogeneous form of the speficied b.c.s

	Assuming (Subscript[ϕ, j]), j=0,N of polynomial form, start with 
	(Subscript[ϕ, 0])(x) = ax +b, from the conditions above, a=1, b=0, (Subscript[ϕ, 0])(x) = x
	(Subscript[ϕ, 1])(x) = x (2-x), which also satisfies the second set of conditions for j>1
	(Subscript[ϕ, 2])(x) = (x^2)(1-(2/3)x), satisfies b.c.s

The approximation function (Subscript[U, 2]) is:

(Subscript[ϕ, 0])= x;
(Subscript[ϕ, 1])= x (2-x);
(Subscript[ϕ, 2])= (x^2)(1-(2/3)x);
(Subscript[U, 2])= (Subscript[ϕ, 0])+(Subscript[c, 1])(Subscript[ϕ, 1])+ (Subscript[c, 2])(Subscript[ϕ, 2])//Simplify

The residual is:

(Subscript[R, 2])= -(∂_(x,x))(Subscript[U, 2]) - (Subscript[U, 2])+ (x^2) //Simplify

Using Galerkin weighting with (Subscript[w, 1])= (Subscript[ϕ, 1])

(Subscript[wr, 1])=(∫_0^1)(Subscript[ϕ, 1]) (Subscript[R, 2])\[DifferentialD]x

Now with (Subscript[w, 2])=(Subscript[ϕ, 2])

(Subscript[wr, 2])=(∫_0^1)(Subscript[ϕ, 2]) (Subscript[R, 2])\[DifferentialD]x

Solving for (Subscript[c, j]),j=1,2

galerkinCs=Solve[{(Subscript[wr, 1])==0, (Subscript[wr, 2])==0},{(Subscript[c, 1]),(Subscript[c, 2])}]//Flatten

(Subscript[UGalerkin, 2])= (Subscript[U, 2])/. galerkinCs//Expand
N[%]

Least squares method
The method consists of finding the set of coefficients (Subscript[c, j]) which minimize the integral of the square of the residual:

		(min_(Subscript[c, j]))((∫_Ω) (([R((Subscript[U, N]))])^2)dΩ)
or
		(∫_Ω) 2 R((Subscript[U, N])) ((∂R((Subscript[U, N])))/(∂(Subscript[c, i])))dΩ = 0, i=1,N
which obviosly corresponds to choosing the weight functions as:

		(Subscript[w, i])=((∂R((Subscript[U, N])))/(∂(Subscript[c, i]))), i=1,N

Example: The two weight functions in this case are

(Subscript[w, 1])= (∂_(Subscript[c, 1]))(Subscript[R, 2])
(Subscript[w, 2])= (∂_(Subscript[c, 2]))(Subscript[R, 2])

The least-squares equations are:

(Subscript[ls, 1])=(∫_0^1) (Subscript[w, 1])(Subscript[R, 2])\[DifferentialD]x
(Subscript[ls, 2])=(∫_0^1) (Subscript[w, 2])(Subscript[R, 2])\[DifferentialD]x

lsCs=Solve[{(Subscript[ls, 1])==0, (Subscript[ls, 2])==0},{(Subscript[c, 1]),(Subscript[c, 2])}]//Flatten

(Subscript[ULSquares, 2])= (Subscript[U, 2])/. lsCs//Expand
N[%]

(Subscript[UGalerkin, 2])//N

The Petrov-Galerkin method: the weight functions are chosen independently of the approximation functions. Going back to the example:

(Subscript[ψ, 1])= x; (Subscript[ψ, 2])= (x^2);

The weighted residual equations are:

(Subscript[pg, 1])=(∫_0^1) (Subscript[ψ, 1])(Subscript[R, 2])\[DifferentialD]x
(Subscript[pg, 2])=(∫_0^1) (Subscript[ψ, 2])(Subscript[R, 2])\[DifferentialD]x

pgCs=Solve[{(Subscript[pg, 1])==0, (Subscript[pg, 2])==0},{(Subscript[c, 1]),(Subscript[c, 2])}]//Flatten

(Subscript[Upg, 2])= (Subscript[U, 2])/. pgCs//Expand
N[%]

(Subscript[UGalerkin, 2])//N

The collocation method: requires the residual to vanish at "N" points inside the domain, e.g. for 0<x<1, N=2, require (Subscript[R, 2])(1/3)=0, (Subscript[R, 2])(2/3) = 0
In the example:

(Subscript[cm, 1])= (Subscript[R, 2])/. x-> 0.25
(Subscript[cm, 2])= (Subscript[R, 2])/. x->0.75
cmCs=Solve[{(Subscript[cm, 1])==0, (Subscript[cm, 2])==0},{(Subscript[c, 1]),(Subscript[c, 2])}]//Flatten

(Subscript[Ucm, 2])= (Subscript[U, 2])/. cmCs//Expand
N[%]

(Subscript[Upg, 2])//N

Now let's get the exact solution, which actually exists for this problem:

exact = DSolve[{-u''[x]-u[x] + (x^2)==0, u[0] ==0, u'[1]==1}, u[x], x]//Flatten;
uu = u[x] /. exact

Plot[{uu, (Subscript[UGalerkin, 2]), (Subscript[ULSquares, 2]),(Subscript[Upg, 2]),(Subscript[Ucm, 2])}, {x,0,1},PlotStyle->{Thick,Red, Green, Blue, Black, Orange},PlotRange-> {{0.9,1},{0.8,1.4}}]

Weak formulation:

Consists of 
1) integrating by parts the weighted residual formulation. 
2) applying the "essential" and "natural" boundary conditions within the formulation

1) Integration by parts
Reminder:
(∫_a^b)v (dw/dx)dx = (∫_a^b)(-(dv/dx)w + ((d(v w))/dx))dx = -(∫_a^b)(dv/dx)wdx + (Subscript[(([v w])^b), a])
This effectively trades derivatives between "u" and "w" and reduces the order of the highest derivative of "u".

Going back to our example:
		(∫_0^1)[w [((-d)/dx)(a (du/dx)))-w f]dx =
(∫_0^1) [(dw/dx)a(du/dx)-w f]dx -(Subscript[(([w a (du/dx)])^1), 0])=
(∫_0^1) [(dw/dx)a(du/dx)-w f]dx -(w a (du/dx))(1) + (w a (du/dx))(0)=0

We say it's a weak formulation of the problem because it now requires  "u" (and (Subscript[U, N])) to be differentiable only once, instead of twice, as in the weighted residual formulation.

2) Application of the boundary conditions:

	- essential boundary conditions: we will require the function "w" to vanish where the primary variable of the problem is specified. In this case:
			u(0) = (Subscript[u, 0]),  then we will require
			w(0) = 0
and the last term in the weak formulation vanishes:
		(∫_0^1) [(dw/dx)a(du/dx)-w f]dx -(w a (du/dx))(1) =0

	- natural boundary conditions: if we apply the boundary condition in the "flux" or secondary variable:
			a(du/dx)(1) = (Subscript[Q, 0]) 
to the weak formulation we obtain:
		(∫_0^1) [(dw/dx)a(du/dx)-w f]dx -w(1) (Subscript[Q, 0]) =0
This is the so-called weak formulation of the problem. It is a restatement of the differential equation with lower continuity requirements on the primary variable AND it also enforces the natural boundary condition.