Green functions with mathematica

vs 0.3, 25/1/07 © Niels Walet

Here we work two problems, and pose a final one, which you should try yourself. Double click the "closed cells" to show the full content. Use the help facility to figure out any unknown constructs.

Problem 1

Find the Green function G(x,x')  for the operator ℒ specified by ℒ y(x)=-((d^2)/(d (x^2)))y(x), in the range 0≤x≤L, with boundary conditions u(0)=u(L)=0. Use as many different techniques as possible.

First we use the eigenvalue method; unfortunately Mathematica needs a little help..

sol=y[x]/.First[DSolve[{-y''[x]==λ y[x],y[0]==0},y[x],x]]

(The "First" gets rid of some curly brackets generated by mathematica, the "/." operator substitutes the result).
Clearly we find

Solve[L (Sqrt[λ])==n π,λ]

sol = a Sin[n π x /L]

Normalise the solutions

(∫_0^L)(sol^2)\[DifferentialD]x==1

Use the fact that n is an integer (which mathematica doesn't know)

Solve[((∫_0^L)(sol^2)\[DifferentialD]x==1/.Sin[2n π]->0),a]

Pick one of these two, and conclude that

G(x,x')=(∑_(n=1)^∞)((L^2)/((n^2) (π^2)))(2/L)sin(((n π x)/L))sin(((n π x')/L))

We use xp rather than x', since Mathematica assigns a special meaning to primes:

Plot3D[(∑_(n=1)^10)((L^2)/(((n )^2)(π^2)))(2/L)Sin[((n π x)/L)]Sin[((n π xp)/L)]/.L->1,{x,0,1},{xp,0,1},AxesLabel->{x,x',G}]

If you so wish this can be simplified; mathematica will even do the sum for you...

Simplify[(∑_(n=1)^∞)((L^2)/(((n )^2)(π^2)))(2/L)Sin[((n π x)/L)]Sin[((n π z)/L)]]

This is a subtle sum of two complex logarithms; if we plot it we see it makes sense, but we have to remove small numerical errors that lead to imaginary parts.

Plot3D[Re[((L ((Subscript[Li, 2])((E^(-((I π (x-z))/L))))+(Subscript[Li, 2])((E^((I π (x-z))/L)))-(Subscript[Li, 2])((E^(-((I π (x+z))/L))))-(Subscript[Li, 2])((E^((I π (x+z))/L)))))/(2 (π^2)))/.L->1],{x,0,1},{z,0,1},AxesLabel->{x,x',G}]

The second method is the continuity one:

First for x<x',

sol=y[x]/.First[DSolve[{-y''[x]==0,y[0]==0},y[x],x]]

Then for x>x',

sol=y[x]/.First[DSolve[{-y''[x]==0,y[L]==0},y[x],x]]

So we have two solutions,`

G(x,x')={({{a(x') x, x<x'}, {b(x') (x-L), x>x'}})

We now use continuity of G, and discontinuity of G' at x=x',

(∫_(x'-ϵ)^(x'+ϵ))-((∂^2)/(∂ (x^2)))G(x,x')dx=1,

or

-(Subscript[[(∂/(∂ x))G(x,x')], x'-ϵ]^(x'+ϵ))=-(∂/(∂ x))G(x'+ϵ,x')+(∂/(∂ x))G(x'-ϵ,x')=1

Don't use x' in Mathematica, use xp instead:

Solve[{a(xp) x==b(xp)(x-L),-(∂_x)(b(xp) (x-L))+(∂_x)(a(xp) x)==1}/.x->xp,{a(xp),b(xp)}]

Plot3D[If[x<xp,((L-xp)/L)x,(xp/L)(L-x)]/.L->1,{x,0,1},{xp,0,1},AxesLabel->{x,x',G}]

Finally we can let Mathematica evaluate the Green function

sol=y[x]/.First[DSolve[{-y''[x]==DiracDelta[x-xp]},y[x],x]]

The θ function or "step function" is zero for negative argument, and 1 for a positive one (another definition of the δ-function: it is the derivative of the θ function)

sol/.x->0

Since θ is 0 for negative argument, (Subscript[c, 1])=0, and thus

sol=x  c-(x-xp) θ(x-xp)

Simplify[sol/.x->L,Assumptions->{xp<L}]

Solve[%==0,c]

Which states c=(L-x')/L,

G(x,y)=(x/L)(L-x')-(x-x')θ(x-x')

Plot3D[((x/L) (L-xp)-(x-xp) UnitStep(x-xp))/.L->1,{x,0,1},{xp,0,1},AxesLabel->{x,x',G}]

Let us check that the answer is correct

First use Mathematica on one problem

DSolve[{-y''[x]==Sin[π x/L],y[0]==0,y[L]==0},y[x],x]

And then integrate the Green function

Simplify[Integrate[(xp/L)(x-L)Sin[π xp/L],{xp,0,x}]+Integrate[-((L-xp)/L)x  Sin[π xp/L],{xp,x,L}]]

Problem 2

Find the Green function G(x,x')  for the operator ℒ y(x)=((d^( ))/(d x))(x (d/dx) y(x)), in the range 0≤x≤1, with the BC that u(1) =0 and u(0) is finite. Use the continuity method and a direct solution from Mathematica.

Continuity method

Gerneral solution with zero r.h.s.

sol=y[x]/.First[DSolve[{(∂_x)(x y'[x])==0},y[x],x]]

(log is ln in Mathematica notation.)

So we have two solutions,

G(x,x')={({{a(x') , x<x' (G(0,x') is finite)}, {b(x') log(x), x>x' (G(1,x')=0)}})

We now use continuity of G, and discontinuity of G' at x=x',

(∫_(x'-ϵ)^(x'+ϵ))(∂/(∂ x))(x (∂/(∂ x))G(x,x'))dx=1, or (Subscript[[x(∂/(∂ x))G(x,x')], x'-ϵ]^(x'+ϵ))=-x'(∂/(∂ x))G(x'+ϵ,x')+x'(∂/(∂ x))G(x'-ϵ,x')=1

Don't use x' in Mathematica, use xp instead:

Solve[{a(xp) ==b(xp)log(x),x(∂_x)(b(xp) log(x))-x(∂_x)(a(xp) )==1}/.x->xp,{a(xp),b(xp)}]

Plot3D[If[x<xp,Log[xp],Log[x]],{x,0.0001,1},{xp,0.0001,1},AxesLabel->{x,x',G}]

Finally we can let Mathematica evaluate the Green function

sol=y[x]/.First[DSolve[{(∂_x)(x y'[x])==DiracDelta[x-xp]},y[x],x]]

sol/.x->0

Since θ is 0 for negative argument, (Subscript[c, 1])=0

sol=c+(log(x)-log(xp)) θ(x-xp)

Solve[Simplify[sol/.x->1,Assumptions->{xp<1}]==0,c]

G(x,x')=log(x')+(log(x)-log(x')) θ(x-x')

Plot3D[log(xp)+(log(x)-log(xp)) θ(x-xp),{x,0.001,1},{xp,0.001,1},AxesLabel->{x,x',G}]

Problem 3

Find the Green function G(x,x')  for the operator ℒ y(x)=((d^( 2))/(d (x^2))) y(x)-(k^2)y(x), in the range -∞≤x≤∞, with boundary conditions y(±∞)=0. Use the eigenvalue and continuity method. Compare to a direct solution from Mathematica.