Re: PDE
- To: mathgroup at smc.vnet.net
- Subject: [mg35346] Re: PDE
- From: Goyder Dr HGD <H.Goyder at rmcs.cranfield.ac.uk>
- Date: Tue, 9 Jul 2002 06:47:59 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
mitja Lakner wrote:
>
> Is there any package for solving eigenvalue problem for the Laplacian
operator
> on the 2D (not rectangular) region with zero boundary value?
>
> f_{xx}+f_{yy}= t f
>
> Mitja Lakner
Mitja,
In my view the inability to solve elliptic PDEs is a major gap in
Mathematica. Below I provide a very poor method of solving these equations
which the experts may be able to improve. However, the method will work for
a non rectangular region which you require.
The eigenvalue problem may be formulated in the standard way. First multiply
both sides of the differential equation by the unknown solution and
integrate over the domain. Green's function applied to the first side
reduces the derivative by one. The problem thus reduces to minimising the
function
-Integrate[Derivative[1, 0][f][#1, #2])^2 + (Derivative[0, 1][f][#1,
#2])^2 & [x,y],dx dy]+
t Integrate[ f[#1,#2]^2, dx dy]
The approach I take below is to approximate the function f with a symbolic
interpolation function. Mathematica is good at doing this. However, although
you can take the derivative of a symbolic interpolation function you cannot
take the integral. Also, to fit a non rectangular region you must put the
domain in a rectangular region and then set all values outside the required
domain to zero.
The worst bit of my formulation below is to approximate the integrals above
by just summing over the domain. I hope someone can advise me on a better
method.
In the example below I do a rectangular region with a square taken out near
the origin.
ClearAll[a]; Lx = 10; Ly = 16; nx = 11; ny = 17; (* define rectangular
region and grid dimensions *)
(* set up zero boundary conditions on edge and outside domain *)
bc = Flatten[Join[
Table[a[i, j] -> 0, {i, 5}, {j, 5}],
Table[a[i, 1] -> 0, {i, 6, nx - 1}],
Table[a[nx, j] -> 0, {j, ny - 1}],
Table[a[i, ny] -> 0, {i, nx - 1, 1, -1}],
Table[a[1, j] -> 0, {j, ny - 1, 6, -1}]]];
(* generate an array of unknown function values in domain *)
f = Table[{Lx(i - 1)/(nx - 1), Ly(j - 1)/(ny - 1), a[i, j]}, {i, nx},
{j,ny}] /. bc;
ff = Interpolation[Flatten[f, 1]]; (* make an interpolation function for
domain *)
d11 = (Derivative[1, 0][ff][#1, #2])^2 + (Derivative[0, 1][ff][#1, #2])^2 &;
(* function to get grad.grad *)
(* estimate the first integral *)
e1 = Apply[Plus,
Flatten[Table[
d11[x, y], {x, 0, Lx, N[Lx/(3nx)]}, {y, 0, Ly, N[Ly/(3ny)]}]]]
// Expand;
(* estimate the second integral *)
e2 = Apply[Plus,
Flatten[Table[
ff[x, y]^2, {x, 0, Lx, N[Lx/(3 nx)]}, {y, 0, Ly,N[Ly/(3 ny)]}]]]
// Expand;
uk = Variables[e1]; (* list of unknowns *)
(* matrices for unknown coefficients *)
m1 = Table[Coefficient[D[e1, uk[[i]]], uk], {i, Length[uk]}];
m2 = Table[Coefficient[D[e2, uk[[i]]], uk], {i, Length[uk]}];
(* standard eigenvalue problem vals are eigenvalues, vecs are eigenvectors
*)
{vals, vecs} = Eigensystem[m1.Inverse[m2]]; vals = Reverse[vals]; vecs =
Reverse[vecs];
ff1 = ff /. Table[uk[[i]] -> vecs[[1, i]], {i, Length[uk]}]; (* make
interpolation function for first eigenfunction *)
ff2 = ff /. Table[uk[[i]] -> vecs[[2, i]], {i, Length[uk]}]; (* make
interpolation function for second eigenfunction *)
Plot3D[ff1[x, y], {x, 0, Lx}, {y, 0, Ly}, PlotPoints -> {50, 50},BoxRatios
-> {Lx, Ly, 5}];
Plot3D[ff2[x, y], {x, 0, Lx}, {y, 0, Ly}, PlotPoints -> {50, 50},BoxRatios
-> {Lx, Ly, 5}];
Not too bad to solve Laplace's equation in 12 lines of code -but the above
could be greatly improved.
Hugh Goyder