L-38 MCS 494 Fri 22 Nov 2002

> restart;

Galerkin's method

We use the example in the book, on pages 244-245.

> ode := diff(y(x),x) + y(x) - 1 = 0;

> init := y(0) = 0; # initial condition

Maple can solve the differential equation symbolically:

> exact := dsolve({ode,init},y(x));

We propose to work with N basis functions x^i, i=1..N.

> N := 5;

The trial solution is then a linear combination of the basis functions:

> trial := sum(c[i]*x^i,i=1..N);

The residual is what is left from the ODE is we plug in the trial solution:

> residual := lhs(subs(y(x)=trial,ode));

> R := simplify(residual);

Galerkin's requirement is that the inner product of the residual with the basis functions is zero. This leads to a linear system in the coefficient of the trial function.

> sys := []:

> for i from 1 to N do

> sys := [op(sys),int(R*x^i,x=0..1) = 0]:

> od:

> sys;

> vars := seq(c[i],i=1..N);

> A,b := LinearAlgebra[GenerateMatrix](sys,[vars]);

> sol := solve({op(sys)},{vars});

> assign(sol);

> trial;

> trial_fun := unapply(trial,x);

> exact_fun := unapply(rhs(exact),x);

> evalf(trial_fun(1) - exact_fun(1));

So the solution is very accurate: on a plot we will not notice the difference between the exact and the approximate solution obtain by Galerkin's method.

The problem with Galerkin's method is that the linear systems become very ill conditioned, i.e.: the solution is highly sensitive to errors on the input, and consequently, the problem is hard to solve with floating-point arithmetic.

> LinearAlgebra[ConditionNumber](A);

> evalf(%);

As the condition number is about 10^6, this means that in order to be sure of having 10 digits correct in the answer, we must use 16 decimal places in our computations.

Note that this is not a problem for Maple, which solved the linear system using exact arithmetic. However, please bear in mind that the solution is a polynomial of degree N. As N increases, the accurate evaluation of the solution will require more and more decimal places.

The Poisson Problem

The following PDE occurs as a model of steady fluid flow down a square tube:

> restart;

> pde := diff(u(x,y),x$2) + diff(u(x,y),y$2) + 1 = 0;

We take zero boundary conditions on the unit square.

As basis functions we consider (x^i-x^(i+1))*(y^j-y^(j+1)).

> N := 4:

> trial := sum(sum(a[i,j]*(x^i-x^(i+1))*(y^j - y^(j+1)),j=1..N),i=1..N);

> residual := lhs(subs(u(x,y)=trial,pde));

> R := simplify(residual);

Imposing Galerkin's requirement leads to N^2 linear equations:

> sys := []:

> for i from 1 to N do

> for j from 1 to N do

> sys := [op(sys),int(int(R*(x^i-x^(i+1))*(y^j-y^(j+1)),x=0..1),y=0..1)]:

> end do:

> end do:

> sys;

> vars := seq(seq(a[i,j],i=1..N),j=1..N);

> sol := solve({op(sys)},{vars});

> assign(sol);

> trial;

How do we test the accuracy of the solution? Well, we evaluate the residual:

> resfun := unapply(simplify(subs(u(x,y)=trial,lhs(pde))),x,y);

> plot3d(resfun,0..1,0..1,axes=boxed);

We see that the solution is not very accurate at the corners. Let us now plot the actual solution:

> plot3d(trial,x=0..1,y=0..1,axes=boxed);

This models the flow of a fluid through a square pipe.

>