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