The Fourier-Laplace Method

For converting an ODE

		p(z) (d^2) (w/d)(z^2) + q(
                  z) d(w/d)z + r(z) w = 
            0

into a relationship

		\(\([F(z, s)]\)\_\(s = a\)\%b + ∫\+C\( e\^\(i\ s\ z\)\) L\ \(f(s)\) = 0\)
wherep(z),q(z),andr(z)are polynomials of arbitrary order.

Description

We seek a solution to the ODE of the form

\(TraditionalForm\`\(\(			\)\(w(
        z)\  = \ ∫\+C\( e\^\(i\ s\ z\)\) \(f(s)\) ds\)\)\)

where C is a contour in the complex plane to be determined.  Begin by computing (w^′)(z) and (w^″)(z).

\(TraditionalForm\`\(\(			\)\(\(w\^′\)(
        z)\  = i\ \(∫\+C 
            s\ \(e\^\(i\ s\ z\)\) \(f(
              s)\) ds\)\)\)\)

\(TraditionalForm\`\(\(			\)\(\(w\^″\)(
        z)\  = \(-\ \(∫\+C s\^2\ \(e\^\(i\ s\ z\)\) \(f(s)\) ds\)\)\)\)\)

Next we want to compute r(z) w(z), q(z) (w^′)(z) and p(z) (w^″)(z).  It will be advantageous to have z appear only in the exponential term of each integrand; other appearances are removed by integrating by parts: \(TraditionalForm\`\(\(∫\_C u\ dv = u\ v\)\(|\)\)\_a\%b - ∫\_C v du\), where a and b denote respectively the initial and terminal points of the contour C.  Later, we will see that one way to achieve our goal is to have a = b, i.e. a closed contour, and so it will not be necessary to keep track of the terms that get evaluated at the endpoints.  Instead, we will denote these terms by ( ... ).

Suppose that   p(
          z)  =  (∑_(j = 0)^l) (Subscript[α, j]) (z^j),    q(
        z)  =  (∑_(j = 0)^m) (Subscript[β, j]) (z^j),   r(
        z)  =  (∑_(j = 0)^n) (Subscript[γ, j]) (z^j).

The following steps show how to compute r(z) w(z).

\(TraditionalForm\`\(\(z\ \(w(z)\)\  = \ \(∫\+C z\ \ \(e\^\(i\ s\ z\)\) \(f(s)\) ds = \((…)\)\)\)\(|\)\)\_\(s = a\)\%b + \(\(\(∫\)\(\ \)\)\+C\) i\ \(e\^\(i\ s\ z\)\) \(d\/ds\) \((\ 
          f(s))\) ds\)

\(TraditionalForm\`\(\(\(\(z\^2\ \(w(z)\)\  = z\ \((…)\)\)\(|\)\)\_\(s = a\)\%b + ∫\+C i\ z\ \(e\^\(i\ s\ z\)\) \(d\/ds\) \((\ f(s))\) ds = \((…)\)\)\(|\)\)\_\(s = a\)\%b - \(\(\(∫\)\(\ \)\)\+C\) \(e\^\(i\ s\ z\)\) \(d\^2\/ds\^2\) \((\ f(s))\) ds\)

\(TraditionalForm\`\(\(\(\(z\^3\ \(\(w\^′\)(z)\)\  = z\ \((…)\)\)\(|\)\)\_\(s = a\)\%b - \(\(\(∫\)\(\ \)\)\+C\) z\ \(e\^\(i\ s\ z\)\) \(d\^2\/ds\^2\) \((\ f(s))\) ds = \((…)\)\)\(|\)\)\_\(s = a\)\%b - ∫\+C i\ \(e\^\(i\ s\ z\)\) \(d\^3\/ds\^3\) \((\ 
            f(s))\) ds\)

\(TraditionalForm\`\(\(z\^k\ \(w(z)\)\  = \((…)\)\)\(|\)\)\_\(s = a\)\%b + \(i\^k\) \(∫\+C\( e\^\(i\ s\ z\)\) \(d\^k\/ds\^k\) \((\ 
              f(s))\) ds\)\)

\(TraditionalForm\`\(\(\(r(z)\)\ \(w(z)\)\  = \((…)\)\)\(|\)\)\_\(s = a\)\%b + ∫\+C\(\( e\^\(i\ s\ z\)\)(∑\+\(k = 0\)\%n\( i\^k\) \(γ\_k\) \(d\^k\/ds\^k\) \((\ 
                  f(s))\))\) ds\)

By similar means we get

\(TraditionalForm\`\(\(\(q(z)\)\ \(\(w\^′\)(z)\)\  = \((…)\)\)\(|\)\)\_\(s = a\)\%b + ∫\+C\(\( e\^\(i\ s\ z\)\)(∑\+\(k = 0\)\%m\( i\^\(k + 
                      1\)\) \(β\_k\) \(d\^k\/ds\^k\) \((s\ \(f(s)\))\))\) ds\)

\(TraditionalForm\`\(\(\(p(z)\)\ \(\(w\^″\)(z)\)\  = \((…)\)\)\(|\)\)\_\(s = a\)\%b + ∫\+C\(\( e\^\(i\ s\ z\)\)(∑\+\(k = 0\)\%l\( i\^\(k + 
                      2\)\) \(α\_k\) \(d\^k\/ds\^k\) \((s\^2\ \(f(s)\))\))\) ds\)

Making the appropriate substitutions into the ODE yields

\(TraditionalForm\`\(\(			\)\((…)\)\(|\)\)\_\(s = a\)\%b + ∫\+C\( e\^\(i\ s\ z\)\) L\ \(f(s)\) ds = 0\)

where L is the linear operator defined by

\(TraditionalForm\`L\ \(f(
          s)\) = ∑\+\(k = 0\)\%l\( i\^\(k + 
                  2\)\) \(α\_k\) \(d\^k\/ds\^k\) \((s\^2\ \(f(
                  s)\))\) + ∑\+\(k = 0\)\%m\( i\^\(k + 
                  1\)\) \(β\_k\) \(d\^k\/ds\^k\) \((s\ \(f(
                  s)\))\) + ∑\+\(k = 0\)\%n\( i\^k\) \(γ\_k\) \(d\^k\/ds\^k\) \((\ 
              f(s))\)\)

One way to achieve \(TraditionalForm\`\(\((…)\)\(|\)\)\_\(s = a\)\%b + ∫\+C\( e\^\(i\ s\ z\)\) 
              L\ \(f(s)\) ds = 0\) is by choosing a = b (i.e. C is a closed contour)  and setting L f(s) = 0.  (For other possible contours see Hille, Ordinary Differential Equations in the Complex Plane [Dover 1997], pp. 216-227.)  That is, we must solve a new ODE which we hope is easier to solve than the original one.  If we succeed, we then choose C by considering the analytic nature of (e^(i s z)) f(
          s).  This notebook is concerned only with arriving at L f(s) and then, if possible, solving for f(s).

Example 1:   z (w^″)(
              z) + (2  ρ + 1) (w^′)(z) + 
          z w(z) = 0

Remove["Global`*"]

Input the left side of the ODE into Mathematica:

Lwz = z w''[z]+(2 ρ+1) w'[z]+z w[z]

Extract the coefficient polynomials:

p[z_]=Coefficient[Lwz,w''[z]];
q[z_]=Coefficient[Lwz,w'[z]];
r[z_]=Coefficient[Lwz,w[z]];
{p[z],q[z],r[z]}

Determine the degrees of the polynomials:

l=Exponent[p[z],z];
m=Exponent[q[z],z];
n=Exponent[r[z],z];
{l,m,n}

Extract the coefficients from each polynomial:

Map[Set[α[#],CoefficientList[p[z],z][[#+1]]]&,Range[0,l]];
Map[Set[β[#],CoefficientList[q[z],z][[#+1]]]&,Range[0,m]];
Map[Set[γ[#],CoefficientList[r[z],z][[#+1]]]&,Range[0,n]];

?α

?β

?γ

Construct L f(s):

sum1 = Sum[I^(k+2) α[k] D[s^2 f[s],{s,k}],{k,0,l}];
sum2 = Sum[I^(k+1) β[k] D[s f[s],{s,k}],{k,0,m}];
sum3 = Sum[I^k γ[k] D[f[s],{s,k}],{k,0,n}];
Lfs = Expand[sum1+sum2+sum3]

Collect terms (optional):

Lfs = Collect[Lfs,Map[D[f[s],{s,#}]&,Range[0,Max[{l,m,n}]]]]

Solve the new ODE:

sol = DSolve[Lfs==0,f[s],s,DSolveConstants->c]

f[s_]=Simplify[f[s]/.sol[[1]]]

Block Form

We can combine the steps from Example 1 into a Mathematica block or module.  This will make it possible to find and solve L f(s) with just one command.  To keep the code simple, the following block requires that the user input the left side of an ODE L w(z) = 0 in terms of w(z) and z.  Having other constants in the coefficient polynomials is allowed.  A good, general-purpose block suitable for inclusion in a Mathematica package would let the user use other variables besides w and z and would also do copious error checking.  Here we just show the rudiments of building a block.  (A module is similar but the local variables are handled differently; see the Mathematica documentation.)  The following block returns a list of two items: the new ODE, and the solution to that ODE.

Remove["Global`*"]

FourierLaplace[Lwz_]:=Block[
	{p,q,r,l,m,n,α,β,γ,sum1,sum2,sum3,Lfs,sol,f,g},
	p[z]=Coefficient[Lwz,w''[z]];
	q[z]=Coefficient[Lwz,w'[z]];
	r[z]=Coefficient[Lwz,w[z]];
	l=Exponent[p[z],z];
	m=Exponent[q[z],z];
	n=Exponent[r[z],z];
	Map[Set[α[#],CoefficientList[p[z],z][[#+1]]]&,Range[0,l]];
	Map[Set[β[#],CoefficientList[q[z],z][[#+1]]]&,Range[0,m]];
	Map[Set[γ[#],CoefficientList[r[z],z][[#+1]]]&,Range[0,n]];
	sum1 = Sum[I^(k+2) α[k] D[s^2 f[s],{s,k}],{k,0,l}];
	sum2 = Sum[I^(k+1) β[k] D[s f[s],{s,k}],{k,0,m}];
	sum3 = Sum[I^k γ[k] D[f[s],{s,k}],{k,0,n}];
	Lfs = Expand[sum1+sum2+sum3];
	Lfs = Collect[Lfs,Map[D[f[s],{s,#}]&,Range[0,Max[{l,m,n}]]]];
	sol = DSolve[Lfs==0,f[s],s,DSolveConstants->c];
	g[s_]=Simplify[f[s]/.sol[[1]]];
	{Lfs==0,g[s]}
	]

Calling our FourierLaplace block with the input from Example 1 yields the expected results:

FourierLaplace[z w''[z]+(2 ρ+1) w'[z]+z w[z]]

Example 2:    (z^2) (w^″)(z) + 
          z (w^′)(z) +  (z^2) w(z) = 0

In Example 1, the new ODE was of first order, i.e. it was simpler than the original.  In this example, we again get a second-order ODE.

FourierLaplace[z^2 w''[z]+ z w'[z]+ z^2 w[z]]

Example 3:   z (w^″)(
              z) + (z^2) (w^′)(z) + (z^3) w(z) = 0

This example shows that it is possible to end up with an ODE of higher order than we started with.

FourierLaplace[z^3 w''[z]+z^2 w'[z]+z w[z]]