Application Center - Maplesoft
  Submit your work
Application Center Applications Green's functions Preview

App Preview:

Green's functions

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application



Introduction to Green's Functions

     Dr. Jack Wagner, Fordham University
     <ajawagner@fordham.edu>                                           copyright 2012

.    Examine the general linear, first order, non-homogenous equation, and its homogeneous counterpart.  

          p(x)*(diff(y(x), x))+g(x)*y(x) = f(x)

          p(x)*(diff(y(x), x))+g(x)*y(x) = 0

    In order to find a solution to the non-homogeneous equation we first look for a solution to the homogeneous one.

With Maple, we calculate.

> 

de:=p(x)*diff(y(x),x)+q(x)*y(x)=0;

p(x)*(diff(y(x), x))+q(x)*y(x) = 0

(1)
> 

h_sol:=dsolve(de,y(x));

y(x) = _C1*exp(Int(-q(x)/p(x), x))

(2)

              Here we have the solution to the first order homogeneous equation. Note that Maple precedes the constant of integration, _C1, with an underscore.  A particular solution to the non-homogeneous equation can be written as a product of  h_sol  and  some other function of x.  We will use this idea to solve the non-homogeneous equation. Let s(x) be the unknown function which multiplies the solution to the homogeneous equation.  We substitute the product s(x)*h_sol for y(x) in the non-homogeneous equation.

> 

De:=lhs(de)=f(x);        

p(x)*(diff(y(x), x))+q(x)*y(x) = f(x)

(3)
> 

De:=subs(y(x)=s(x)*rhs(h_sol),De);

p(x)*(diff(s(x)*_C1*exp(Int(-q(x)/p(x), x)), x))+q(x)*s(x)*_C1*exp(Int(-q(x)/p(x), x)) = f(x)

(4)

     Now solve for s(x).

> 

s_sol:=dsolve(De,s(x));

s(x) = Int(f(x)*exp(Int(q(x)/p(x), x))/(_C1*p(x)), x)+_C2

(5)

     This is the function which multiplies the solution to the homogeneous equation to produce a particular solution of the non-homogeneous equation.  Because we have already evaluated the homogeneous equation solution, we replace its variable with a "dummy"  variable in order to avoid confusion with the variable with respect to which we are yet to integrate.

> 

h_sol:=subs(x=s,h_sol);

y(s) = _C1*exp(Int(-q(s)/p(s), s))

(6)

      Now let's have a look at our  solution

> 

expand(rhs(h_sol)*rhs(s_sol));

(Int(f(x)*exp(Int(q(x)/p(x), x))/p(x), x))/exp(Int(q(s)/p(s), s))+_C1*_C2/exp(Int(q(s)/p(s), s))

(7)

       What previously appeared as two separate constants of integration have now merged into one as they should.  Now recollect the solution to the homogeneous equation, exp(-(Int(q(s)/p(s), s))) , and define:  

    G[s, x] = y(s)/(y(x)*p(x))

> 

G[s,x]:=exp(-Int(q(s)/p(s),s))/(exp(-Int(q(x)/p(x),x))*p(x));

exp(-(Int(q(s)/p(s), s)))/(exp(-(Int(q(x)/p(x), x)))*p(x))

(8)

       Note that a particular solution of the non-homogeneous equation may be written as:

Int(exp(-(Int(q(s)/p(s), s)))*f(x)/(exp(-(Int(q(x)/p(x), x)))*p(x)), x) =  Int(y(s)*f(x)/(y(x)*p(x)), x) = Int(G[s, x]*f(x), x) 

         

       A differential equation such as the general linear first order differential equation may be regarded as the operation of a linear operator, L= p(x)*d/dx+q(x)  on the solution function,  y(x).  L*y(x) = f(x). The Green's function G[s, x] is then a kind of inverse of L.  y(x)= Int(G[s, x]*L*y(x), x) = Int(G[s, x]*f(x), x)

     In Maple it looks like this.

> 

part_sol:=Int(G[s,x]*f(x),x);

Int(exp(-(Int(q(s)/p(s), s)))*f(x)/(exp(-(Int(q(x)/p(x), x)))*p(x)), x)

(9)

       Once we know the Green's function for a differential operator L,  we may solve the equation L*y(x) = f(x) for any function on the right side by substitution in the above expression.  

Ex. 1

     Find the general solution of  sin(x)*(diff(y(x), x))+cos(x)*y(x) = 1

> 

restart:

> 

De:=sin(x)*diff(y(x),x)+cos(x)*y(x)=1;

sin(x)*(diff(y(x), x))+cos(x)*y(x) = 1

(10)

     We will now apply the Green's function machinery to this equation. First, recollect the definition of the Green's function for first order equations.

> 

G[s,x]:=exp(-Int(q(s)/p(s),s))/(exp(-Int(q(x)/p(x),x))*p(x));

exp(-(Int(q(s)/p(s), s)))/(exp(-(Int(q(x)/p(x), x)))*p(x))

(11)

      Now define p and q for our particular problem.

> 

PQ:={p(x)=sin(x),p(s)=sin(s),q(x)=cos(x),q(s)=cos(s)};

{p(s) = sin(s), p(x) = sin(x), q(s) = cos(s), q(x) = cos(x)}

(12)
> 

G[s,x]:=subs(PQ,G[s,x]);

exp(-(Int(cos(s)/sin(s), s)))/(exp(-(Int(cos(x)/sin(x), x)))*sin(x))

(13)

    Find Int(G[s, x]*f(x), x)     

> 

p_sol:=Int(G[s,x]*1,x);

Int(exp(-(Int(cos(s)/sin(s), s)))/(exp(-(Int(cos(x)/sin(x), x)))*sin(x)), x)

(14)
> 

p_sol:=value(%);

x/sin(s)

(15)

     Replace the "dummy" variable with the independent variable.

> 

p_sol:=subs(s=x,%);

x/sin(x)

(16)

     Recall the solution to the homogeneous equation.

> 

h_sol := _C1*exp(Int(-q(x)/p(x),x));

_C1*exp(Int(-q(x)/p(x), x))

(17)
> 

h_sol:=subs(PQ,h_sol);

_C1*exp(Int(-cos(x)/sin(x), x))

(18)
> 

h_sol:=value(%);

_C1/sin(x)

(19)

The general solution is the sum of the homogeneous and any particular solution.

> 

g_sol:=simplify(h_sol+p_sol);

(_C1+x)/sin(x)

(20)

      If one were to do this with any frequency it would pay to define a procedure that creates the appropriate Green's function.

> 

restart:

> 

G:=proc(p::algebraic,q::algebraic)
    local px,qx,ps,qs,Gxs:              qx:=q: px:=p: ps:=subs(x=s,p): qs:=subs(x=s,q):
    Gxs:=exp(-Int(qs/ps,s))/(exp(-Int(qx/px,x))*px):
end proc:

    

Ex. 2

  We can then use the previously defined procedure (assuming it has been appropriately saved).

> 

de:=x^2*diff(y(x),x)+x*y(x)=sin(x);

x^2*(diff(y(x), x))+x*y(x) = sin(x)

(21)
> 

G:=G(x^2,x);

exp(-(Int(1/s, s)))/(exp(-(Int(1/x, x)))*x^2)

(22)
> 

p_sol:=int(G*sin(x),x);

int(exp(-(Int(1/s, s)))*sin(x)/(exp(-(Int(1/x, x)))*x^2), x)

(23)
> 

p_sol:=value(%);

Si(x)/s

(24)
> 

p_sol:=subs(s=x,p_sol);

Si(x)/x

(25)

    The solution to the homogeneous equation is conveniently found in the numerator of the Green's function.

> 

h_sol:=k*numer(G);

k*exp(-(Int(1/s, s)))

(26)
> 

h_sol:=value(%);

k/s

(27)
> 

h_sol:=subs(s=x,h_sol);

k/x

(28)
> 

g_sol:=simplify(h_sol+p_sol);

(k+Si(x))/x

(29)

       We have gone to an immense amount of trouble to end up with the same answer that we could have obtained  with a very simple Maple function.  

> 

dsolve(de,y(x));

y(x) = (Si(x)+_C1)/x

(30)

      The point, of course, is not the answer but the  introduction to Green's functions and the manner in which they are used.  The importance of this method and the theory that surrounds it can hardly be overstated.

       For something a bit more complcated consider the second order linear differential equation and its homogeneous counterpart.

a[2](t)*(Diff(y(t), t, t))+a[1](t)*(Diff(y(t), t))+a[0](t)*y(t) = f(t)

a[2](t)*(Diff(y(t), t, t))+a[1](t)*(Diff(y(t), t))+a[0](t)*y(t) = 0

      To find the general solution we again start with the homogeneous equation.

> 

restart:

> 

eq1:=a[2](t)*Diff(y(t),t$2)+a[1](t)*Diff(y(t),t)+a[0](t)*y(t)=0;

a[2](t)*(Diff(y(t), t, t))+a[1](t)*(Diff(y(t), t))+a[0](t)*y(t) = 0

(31)

     Now we need solutions to the homogeneous equation, y[1](t) and  y[2](t).

     A particular solution of the non-homogeneous equation can be expressed as:

y[p](t) = u[1](t)*y[1](t)+u[2](t)*y[2](t) for independent functions, u[1](t), u[2](t).

> 

eq2:=y[p](t)=u[1](t)*y[1](t)+u[2](t)*y[2](t);

y[p](t) = u[1](t)*y[1](t)+u[2](t)*y[2](t)

(32)

    Compute the first derivative of  y[p](t).

> 

eq3:=map(diff,eq2,t);

diff(y[p](t), t) = (diff(u[1](t), t))*y[1](t)+u[1](t)*(diff(y[1](t), t))+(diff(u[2](t), t))*y[2](t)+u[2](t)*(diff(y[2](t), t))

(33)

     In addition to the conditions on u[1](t) and u[2](t) expressed in eq2, we may place a second condition, i.e.

> 

eq4:=0=diff(u[1](t),t)*y[1](t)+diff(u[2](t),t)*y[2](t);

0 = (diff(u[1](t), t))*y[1](t)+(diff(u[2](t), t))*y[2](t)

(34)

     This second condition yields.

> 

eq5:=eq3-eq4;

diff(y[p](t), t) = u[1](t)*(diff(y[1](t), t))+u[2](t)*(diff(y[2](t), t))

(35)

     Substitute this expression for diff(y[p](t), t) into the general non-homogeneous second order equation.

> 

eq6:=a[2](t)*Diff(rhs(eq5),t)+a[1](t)*rhs(eq5)+a[0](t)*rhs(eq2)=f(t);

a[2](t)*(Diff(u[1](t)*(diff(y[1](t), t))+u[2](t)*(diff(y[2](t), t)), t))+a[1](t)*(u[1](t)*(diff(y[1](t), t))+u[2](t)*(diff(y[2](t), t)))+a[0](t)*(u[1](t)*y[1](t)+u[2](t)*y[2](t)) = f(t)

(36)
> 

eq7:=value(eq6);

a[2](t)*((diff(u[1](t), t))*(diff(y[1](t), t))+u[1](t)*(diff(diff(y[1](t), t), t))+(diff(u[2](t), t))*(diff(y[2](t), t))+u[2](t)*(diff(diff(y[2](t), t), t)))+a[1](t)*(u[1](t)*(diff(y[1](t), t))+u[2](t)*(diff(y[2](t), t)))+a[0](t)*(u[1](t)*y[1](t)+u[2](t)*y[2](t)) = f(t)

(37)

     Because  y[1](t) and   y[2](t)  are both solutions of the homogeneous equation, we may write:

> 

eq8a:=u[1](t)*subs(y(t)=y[1](t),eq1);

u[1](t)*(a[2](t)*(Diff(y[1](t), t, t))+a[1](t)*(Diff(y[1](t), t))+a[0](t)*y[1](t)) = 0

(38)
> 

eq8b:=u[2](t)*subs(y(t)=y[2](t),eq1);

u[2](t)*(a[2](t)*(Diff(y[2](t), t, t))+a[1](t)*(Diff(y[2](t), t))+a[0](t)*y[2](t)) = 0

(39)
> 

eq9:=simplify(eq7-eq8a-eq8b);

a[2](t)*((diff(u[1](t), t))*(diff(y[1](t), t))+(diff(u[2](t), t))*(diff(y[2](t), t))) = f(t)

(40)

     Recall the second condition on diff(u[1](t), t) and diff(u[2](t), t).

> 

eq4:=eq4;

0 = (diff(u[1](t), t))*y[1](t)+(diff(u[2](t), t))*y[2](t)

(41)

     solve eq9 and eq4 for diff(u[1](t), t) and diff(u[2](t), t).

> 

sol:=solve({eq9,eq4},{diff(u[1](t),t),diff(u[2](t),t)});

{diff(u[1](t), t) = f(t)*y[2](t)/(a[2](t)*(y[2](t)*(diff(y[1](t), t))-(diff(y[2](t), t))*y[1](t))), diff(u[2](t), t) = -y[1](t)*f(t)/(a[2](t)*(y[2](t)*(diff(y[1](t), t))-(diff(y[2](t), t))*y[1](t)))}

(42)

     Now observe that in the denominator of each solution there is the determinant of the Wronskian.

Wronskian

> 

with(VectorCalculus):with(LinearAlgebra):

> 

Wronskian([y[1](t),y[2](t)],t);

Matrix(%id = 18446744078110741798)

(43)
> 

Determinant(%);

(diff(y[2](t), t))*y[1](t)-y[2](t)*(diff(y[1](t), t))

(44)

    As before we are going to need a "dummy" variable.

> 

U[1]:=subs(t=s,Int(rhs(sol[1]),t));

Int(f(s)*y[2](s)/(a[2](s)*(y[2](s)*(diff(y[1](s), s))-(diff(y[2](s), s))*y[1](s))), s)

(45)
> 

U[2]:=subs(t=s,Int(rhs(sol[2]),t));

Int(-y[1](s)*f(s)/(a[2](s)*(y[2](s)*(diff(y[1](s), s))-(diff(y[2](s), s))*y[1](s))), s)

(46)

     

> 

y[p](t):=subs({u[1](t)=U[1],u[2](t)=U[2]},rhs(eq2));

(Int(f(s)*y[2](s)/(a[2](s)*(y[2](s)*(diff(y[1](s), s))-(diff(y[2](s), s))*y[1](s))), s))*y[1](t)+(Int(-y[1](s)*f(s)/(a[2](s)*(y[2](s)*(diff(y[1](s), s))-(diff(y[2](s), s))*y[1](s))), s))*y[2](t)

(47)
> 

y[p](t):=combine(y[p](t));

Int(-y[2](t)*y[1](s)*f(s)/(a[2](s)*(y[2](s)*(diff(y[1](s), s))-(diff(y[2](s), s))*y[1](s)))+y[1](t)*f(s)*y[2](s)/(a[2](s)*(y[2](s)*(diff(y[1](s), s))-(diff(y[2](s), s))*y[1](s))), s)

(48)

      Now define the Green's function.

> 

G[s,t]:=-y[1](t)*y[2](s)/a[2](s)/(diff(y[2](s),s)*y[1](s)-y[2](s)*diff(y[1](s),s))+y[2](t)*y[1](s)/a[1](s)/(diff(y[2](s),s)*y[1](s)-y[2](s)*diff(y[1](s),s));

-y[1](t)*y[2](s)/(a[2](s)*((diff(y[2](s), s))*y[1](s)-y[2](s)*(diff(y[1](s), s))))+y[2](t)*y[1](s)/(a[1](s)*((diff(y[2](s), s))*y[1](s)-y[2](s)*(diff(y[1](s), s))))

(49)

     Having found the Green's function we may solve the non-homogeneous equation for any function on the right hand side.

> 

y[p](t):=Int(expand(G[s,t]*f(s)),s);

Int(-f(s)*y[1](t)*y[2](s)/(a[2](s)*((diff(y[2](s), s))*y[1](s)-y[2](s)*(diff(y[1](s), s))))+f(s)*y[2](t)*y[1](s)/(a[1](s)*((diff(y[2](s), s))*y[1](s)-y[2](s)*(diff(y[1](s), s)))), s)

(50)
> 

y(t):=y[p](t)+c[1]*y[1](t)+c[2]*y[2](t);

Int(-f(s)*y[1](t)*y[2](s)/(a[2](s)*((diff(y[2](s), s))*y[1](s)-y[2](s)*(diff(y[1](s), s))))+f(s)*y[2](t)*y[1](s)/(a[1](s)*((diff(y[2](s), s))*y[1](s)-y[2](s)*(diff(y[1](s), s)))), s)+c[1]*y[1](t)+c[2]*y[2](t)

(51)

 

Ex. 3

    Find the general solution to:         t^2*(Diff(y(t), t, t))+t*(Diff(y(t), t))-y(t) = t

> 

restart: with(VectorCalculus): with(LinearAlgebra):

> 

eq1:=t^2*Diff(y(t),`$`(t,2))+t*Diff(y(t),t)-y(t) = t;

t^2*(Diff(y(t), t, t))+t*(Diff(y(t), t))-y(t) = t

(52)
> 

eq2:=t^2*Diff(y(t),`$`(t,2))+t*Diff(y(t),t)-y(t) =0;

t^2*(Diff(y(t), t, t))+t*(Diff(y(t), t))-y(t) = 0

(53)

      The solution will be a polynomial of low degree.

> 

subs(y(t)=_C1*t+_C2*t^2+_C3*t^3,eq2);

t^2*(Diff(_C1*t+_C2*t^2+_C3*t^3, t, t))+t*(Diff(_C1*t+_C2*t^2+_C3*t^3, t))-_C1*t-_C2*t^2-_C3*t^3 = 0

(54)
> 

expand(value(%));

3*_C2*t^2+8*_C3*t^3 = 0

(55)

    _C3 and C2 must be zero and _C1 is indeterminate.

> 

subs(y(t)=_C4/t+_C5/t^2+_C6/t^3,eq2);

t^2*(Diff(_C4/t+_C5/t^2+_C6/t^3, t, t))+t*(Diff(_C4/t+_C5/t^2+_C6/t^3, t))-_C4/t-_C5/t^2-_C6/t^3 = 0

(56)
> 

expand(value(%));

3*_C5/t^2+8*_C6/t^3 = 0

(57)

_C5 and _C6 must be zero and _C4 is indeterminate.

> 

y[1](t):=_C1*t; y[2](t):=_C4/t;

_C1*t

_C4/t

(58)
> 

y[1](s):=subs(t=s,y[1](t)); y[2](s):=subs(t=s,y[2](t));

_C1*s

_C4/s

(59)
> 

a[2]:=s^2:

macro

> 

macro(W=Determinant(Wronskian([y[1](s),y[2](s)],s))):

> 

W;

-2*_C1*_C4/s

(60)

     Define G[s, t].

> 

G[s,t]:=-y[1](t)*y[2](s)/(a[2]*W)+y[2](t)*y[1](s)/(a[2]*W);

(1/2)*t/s^2-(1/2)/t

(61)

     The particular solution is:

> 

y[p](t):=Int(G[s,t]*s,s);

Int(((1/2)*t/s^2-(1/2)/t)*s, s)

(62)
> 

y[p](t):=value(%);

-(1/4)*s^2/t+(1/2)*t*ln(s)

(63)

     Check the solution.

> 

y[p](t):=subs(s=t,y[p](t));

-(1/4)*t+(1/2)*t*ln(t)

(64)
> 

subs(y(t)=y[p](t),eq1);

t^2*(Diff(-(1/4)*t+(1/2)*t*ln(t), t, t))+t*(Diff(-(1/4)*t+(1/2)*t*ln(t), t))+(1/4)*t-(1/2)*t*ln(t) = t

(65)
> 

value(%);

(3/4)*t+t*(1/4+(1/2)*ln(t))-(1/2)*t*ln(t) = t

(66)
> 

simplify(%);

t = t

(67)
> 

gen_sol:=y[p](t)+y[1](t)+y[2](t);

-(1/4)*t+(1/2)*t*ln(t)+_C1*t+_C4/t

(68)

Ex. 4

     Now it is a simple matter to obtain solutions for a whole family of  equations with varying "impulse" terms. Let k vary from 1 to 5.

      t^2*(Diff(y(t), t, t))+t*(Diff(y(t), t))-y(t) = cos(t)^k

> 

for k from 1 to 5 do
y||k:=Int(G[s,t]*cos(s)^k,s):
y||k:=value(%):od:
for k from 1 to 5 do
y||k:=subs(s=t,y||k); od;

(1/2)*t*(-cos(t)/t-Si(t))-(1/2)*sin(t)/t

(1/2)*t*(-(1/2)*cos(2*t)/t-Si(2*t)-(1/2)/t)-(1/2)*((1/2)*cos(t)*sin(t)+(1/2)*t)/t

(1/2)*t*(-(1/4)*cos(3*t)/t-(3/4)*Si(3*t)-(3/4)*cos(t)/t-(3/4)*Si(t))-(1/2)*((1/3)*cos(t)^2*sin(t)+(2/3)*sin(t))/t

(1/2)*t*(-(1/8)*cos(4*t)/t-(1/2)*Si(4*t)-(1/2)*cos(2*t)/t-Si(2*t)-(3/8)/t)-(1/2)*((1/4)*cos(t)^3*sin(t)+(3/8)*cos(t)*sin(t)+(3/8)*t)/t

(1/2)*t*(-(1/16)*cos(5*t)/t-(5/16)*Si(5*t)-(5/16)*cos(3*t)/t-(15/16)*Si(3*t)-(5/8)*cos(t)/t-(5/8)*Si(t))-(1/2)*((1/5)*cos(t)^4*sin(t)+(4/15)*cos(t)^2*sin(t)+(8/15)*sin(t))/t

(69)

     For convenience in making a comparison we will set _C1 and _C4 equal to one.

> 

g:=t+1/t;

t+1/t

(70)

     The general solution to the original equation is then, g+yk.

> 

plot([g+y1,g+y2,g+y3,g+y4,g+y5],t=0..5,y=-3..10,color=[red,blue,green,black,gold]);

Ex. 5

      t^2*(Diff(y(t), t, t))+t*(Diff(y(t), t))-y(t) = k*exp(k*t)

> 

for k from 1 to 5 do
y||k:=Int(G[s,t]*k*exp(k*s),s):
y||k:=value(%):od:
for k from 1 to 5 do
y||k:=subs(s=t,y||k); od;

(1/2)*t*(-exp(t)/t-Ei(1, -t))-(1/2)*exp(t)/t

2*t*(-(1/2)*exp(2*t)/t-Ei(1, -2*t))-(1/2)*exp(2*t)/t

(9/2)*t*(-(1/3)*exp(3*t)/t-Ei(1, -3*t))-(1/2)*exp(3*t)/t

8*t*(-(1/4)*exp(4*t)/t-Ei(1, -4*t))-(1/2)*exp(4*t)/t

(25/2)*t*(-(1/5)*exp(5*t)/t-Ei(1, -5*t))-(1/2)*exp(5*t)/t

(71)

     Check a typical solution.

> 

eqa:=subs(y(t)=y5,lhs(eq1));

t^2*(Diff((25/2)*t*(-(1/5)*exp(5*t)/t-Ei(1, -5*t))-(1/2)*exp(5*t)/t, t, t))+t*(Diff((25/2)*t*(-(1/5)*exp(5*t)/t-Ei(1, -5*t))-(1/2)*exp(5*t)/t, t))-(25/2)*t*(-(1/5)*exp(5*t)/t-Ei(1, -5*t))+(1/2)*exp(5*t)/t

(72)
> 

value(%);

t^2*(-exp(5*t)/t^3+5*exp(5*t)/t^2)+t*(-(25/2)*Ei(1, -5*t)+(1/2)*exp(5*t)/t^2-(5/2)*exp(5*t)/t)-(25/2)*t*(-(1/5)*exp(5*t)/t-Ei(1, -5*t))+(1/2)*exp(5*t)/t

(73)
> 

simplify(%);

5*exp(5*t)

(74)

Ei

complexplot

> 

plots[complexplot]([g+y1,g+y2,g+y3,g+y4,g+y5],t=0..0.1,view=[0..1000,0..1]);

> 

%?

> 

%?