Web Page Started November 20, 1995 Latest Update November 14, 2019

**INTRODUCTION:
Over the past thirty five years I have taught an
introductory sequence of graduate level courses in
mathematical analysis for engineering majors. An estimated 3000
students have attended these courses and the topics covered
are Advanced Ordinary Differential
Equations (EGM6321), Partial Differential Equations and
Boundary Value Problems(EGM6222), and Integral Equations and
Variational Methods(EGM6323). We
present here a collection of equations and graphs for
the various mathematical functions encountered in our lectures
and give some relevant historical information on the
mathematicians associated with their discovery. Equations
and figures are presented in sequential order as they arise in
the lectures. You can view the
graphs by clicking on the underlined titles or the activated
word HERE. The graphs have been generated by a variety
of canned programs including Maple, Matlab, Mathematica, and
Mathcad. OUTLINES for the three
courses are found HERE,
HERE,andHERE.
In the descriptions below we are using a type of mathematical
esperanto since html does not have the standard math
symbols in its repertoire.**

**You can contact me at E-Mail kurzweg@ufl.edu
or reach me via
standard mail at U.H.Kurzweg, MAE-A Bldg, University of
Florida, Gainesville, FL 32611, USA. This page is found on the internet
at http://www2.mae.ufl.edu/~uhk/MATHFUNC.htm**

**In case you are wondering what an image
of Duerer's etching "Melancholia" is doing on our title
page? It's mainly due to the 4x4 magic square shown above the
head of the seated thinker. Known as Duerer's magic square ,
its matrix reads [16,3,2,13;5,10,11,8;9,6,7,12;4,15, 14,1]
with the semicolons separating rows. Note that the
juxtaposition of elements a _{4,2}=15 and a_{4,3}=14
in the fourth row give the 1514 completion date of the
etching. The numbers along any vertical, horizontal or
diagonal line in his magic square add up to n(n^2+1)/2=34,
which when read backwards gives the artist's age at the time
of the etching. Duerer the artist (1471-1528) obviously also
had more than just a cursory interest in mathematics and
numerology!**

**At
the
bottom of this WEB page you will also find some of our latest thoughts and discussions
on additional mathematical subjects of common
interest not covered directly by our three analysis courses.
They appear in chronological order with the latest at the
end. Some of these topics might be of interest to the
readers as they contain a lot of new and original
observations of possible use for
future studies especially in the area of number theory.
**

**EGM 6321-ADVANCED
ORDINARY DIFFERENTIAL EQUATIONS**

*Solution of linear and non-linear ordinary differential
equations. Method of Frobenius, classification of singularities.*

*Integral representations of solutions. Treatment of the Bessel,
Hermite, Legendre, Hypergeometric and Mathieu*

*equations. Asymptotic methods including the WBK and saddle
point techniques. Treatment of non-linear autonomous*

*equations. Phase plane trajectories and limit cycles.
Thomas-Fermi, Emden, and van der Pol equations.*

**ISOCLINE
SOLUTION OF FIRST ORDER ODE'S: **-** There is available on the WEB a numerical
program written by John Polking for MATLAB ( http://math.rice.edu/~dfield/) which can draw isoclines for the solution of
any first order ODE and then allow one to graph the solution
curve which satisfies a particular initial condition. I
demonstrate the use of this program for the deceptively simple
equation y'=1/x-1/y which clearly has solutions with infinite
slope along both the y and x axis. Can anyone find an analytic
solution to this equation? I've been waiting for one for ten
years.**

**THE LOGISTIC EQUATION: -One of the better known first order ODE's is
the logistic equation y'=y(1-y) which arises in describing
population growth of a species in an environment with limited
resources and also arises in connection with certain studies
in chaos. It is a special case of the Bernoulli equation
y'+p(x)*y=q(x)*y^n and has a straight forward separation of
variables solution y=1/[1+((1-A)/A))exp(-x)] , where y(0)=A
and y(infinity)=1. We give here the graphs of this solution
and its isoclines for several different initial conditions A.**

**OBLIQUE
TRAJECTORIES:- Consider a
family of curves defined by the first order differential
equation y'=f(x,y). Now think of a second family of curves
which intersect these curves at angle a.
These later curves are referred to as the oblique trajectories
and they have the slope y'=(f+m)/(1-f*m) as is readily
established by use of the double angle formula for tangent.
The constant m=tan(a). Consider now
the case of the oblique trajectories intersecting the
concentric circles x^2+y^2=C^2 at 45 degrees. Here f=-x/y and
m=-1, so that the governing equation for the trajectories
becomes y'=(x+y)/(x-y). This is a homogeneous equation
solvable by the variable substitution v=y/x and yields
the logarithmic spiral which (in polar coordinates) reads
ln(r)=[A+B*q]. A plot of one of
these spirals for A=0 and B=1 is shown in the accompanying
graph. The famous mathematician Jacob Bernoulli (1654-1705)
was so enamored with this logarithmic spiral and its
properties that he had it engraved on his tombstone in Basel,
Switzerland. By clicking HERE you can see me pointing to this spiral as it
appears on his epitaph. Sorry about the quality of the jpg. In
the same church you can also find the much more visited grave
of Erasmus of Rotterdam, the famous Dutch theologian and
humanist of the early sixteenth century.**

**TWELVE
IMPORTANT SECOND ORDER LINEAR ODE'S: -We show here a list of the 12 most
important classical second order variable coefficient
linear equations with which most physicists and engineers
should be familiar by the time they reach the masters degree
level. The equations are encountered in a variety of different
areas including optics, quantum mechanics, laser physics,
electromagnetic wave propagation, solid and fluid mechanics
and heat transfer. We will give detailed solutions to each in
the present course and will derive some of the more
interesting generating functions and integral representations
for the solutions .**

**AIRY
FUNCTIONS OF THE FIRST [Ai(x)] AND SECOND [Bi(x)] KIND:
-These follow from a series solution to y"-xy=0.
Note their exponential behaviour for x>0 and oscillatory
character for x<0. They can also be written as Bessel
functions of 1/3 order. The functions are named after George
Airy who was a professor at Cambridge and the Astronomer
Royal of England from 1835 until 1881. Click HERE to see an image of Airy. He
was a bit of a "stuffed shirt" in his position as Astronomer
Royal and failed to follow up with observations on the location
of Neptune predicted by Adams (thus giving LeVerrier and hence
France, and not England, priority for the discovery of this
planet) and also made the comment "of what possible use could
such a machine have" when answering a government inquiry about
supporting the work on the analytical engine(ie.-computer) by
Charles Babbage. On the more positive side, his name is also
attached to the Airy stress function in elasticity and he was
responsible for establishing the Airy prime meridian at
Greenwich. The two convergent infinite series which satisfy the
Airy equation are f(x)=1+x^3/(2*3)+x^6/(2*3*5*6)+ and
g(x)=x+x^4/(3*4)+x^7/(3*4*6*7)+ and the Airy functions are
defined as the linear combinations Ai(x)=a*f(x)-b*g(x) and
Bi(x)=sqrt(3)[a*f(x)+b*g(x)], where a=Ai(0)=0.355028 and
b=-Ai(0)'=0.258819. The functions are tabulated in the Handbook of Mathematical Functions
by Abramowitz and Stegun(Dover Pub) and are contained in canned
programs such as MAPLE and MATHEMATICA.**

**COMPARISON
OF AN ANALYTIC AND A NUMERICAL SOLUTION OF Y"-X*Y=0
SUBJECTED TO Y(0)=0,Y'(0)=1: In view of today's class discussion we know
that y"-x*y=0 has an analytic solution y(x)=a*Ai(x)+b*Bi(x).
When we impose the conditions y(0)=0 and y'(0)=1, the
constants a and b can be evaluated to yield the analytic
result-**

**
y(x)=[-Bi(0)*Ai(x)+Ai(0)*Bi(x)]/[Ai(0)*Bi'(0)-Ai'(0)*Bi(0)] .**

**We have plotted this solution over the
range -8<x<3 in red. To see how this agrees with a
numerical solution using the Runga-Kutta method one can go to
MAPLE . Here the one liner-**

**DEplot(diff(y(x),x$2)=x*y(x),y(x),x=-6..2,[[y(0)=0,D(y)(0)=1]],y=-5..20,linecolour=blue,stepsize=0.02)
;**

**produces the blue curve in the range
-6<x<2 shown. As expected the two solutions yield
identical results. The advantage of an analytic approach(when
this is possible) is that one can obtain a general solution
not dependent on initial or boundary conditions.**

**HERMITE
POLYNOMIALS: -We have shown via a series expansion about
x=0 that the Hermite equation y"-2xy'+2ny=0 has a polynomial
solution whenever n is an integer. The second linearly
independent solution remains an infinite series. These so-called
Hermite polynomials can also be generated (as we will show later
in the course) by the generating function
H[n+2,x]=2xH[n+1,x]-2(n+1)H[n,x], where the term in the square
bracket means a subscript. Starting with H[0,x]=1 and H[1,x]=2x
one finds H[2,x]=4x^2-2, H[3,x]=8x^3-12x , etc. Click HERE
to see a few more of the H(n,x)s computer generated by this
formula. Note the the Hermite polynomials are even when n is
even and odd when n is odd. The orthogonality relation for
H[n,x] polynomials is Int[exp(-x^2)H[n,x]H[m,x],
x=-infinity..+infinity]=2^n n! sqrt(p)delk,
wheredelk
is the Kronecker delta equal to 1 when n=m and zero when integer
n does not equal m. Charles
Hermite(1822-1901) was a mathematician working at the
Ecole Polytechnique and the Sorbonne and well known for his
differential equation, their polynomial solutions, Hermitian
matrixes, work on the quintic equation and its relation to
elliptic functions, plus the proof that e is a transcendental
number.**

**WAVE
FUNCTIONS FOR THE HARMONIC OSCILLATOR: -In class we developed the series solutions
for the Hermite equation y"-2xy'+2n*y=0 and showed that
whenver n is an integer then one of its two solutions is a
Hermite polynomial H(n,x)=(-1)^n*exp(x^2)*D[exp(-x^2), n]. An
important application of these polynomials arises in
connection with solving the Schroedinger equation y"+[2m/(hbar^2)]*[E-V(x)]*y=0 for a harmonic oscillator where the
potential goes as V(x)=(1/2)*k*x^2. As shown in class, the
solution , satisfying the Schroedinger equation and one which
vanishes at plus and minus infinity, is y=H(n,X)*exp-(0.5*X^2), where
X=const*x. A plot for the first three solutions corresponding
to n=0, 1 and 2 is shown. The corresponding eigenvalues are
found to be E=sqrt(k/m)*hbar*(2n+1)/2. The energy levels E are
thus seen to be quantized with a zero point energy of
(1/2)hbar*w=(1/2) h*n. The value of Plank's constant is
h=hbar*(2*p)=6.6256*10^(-34) J sec
and n =w/(2*p) is the oscillator frequency
expressed in Hz.**

**BESSEL
FUNCTIONS OF THE FIRST KIND[J(n,x)]: -The first three Bessel functions of
the first kind corresponding to n=0, 1, and 2. They represent
solutions to x^2y"+xy'+(x^2-n^2)y=0 and are encountered when
formulating certain physical problems in terms of cylindrical
coordinates. The series form for the Bessel functions
J(n,x) is
Sum[(-1)^k*(x/2)^(2k+n)/[k!(n+k)!],{k,0,Inf}]. Friedrich
Bessel was director of the astronomical observatory in
Koenigsberg, East Prussia(now Kaliningrad, Russia) during the
first half of the 19th century. He made the first stellar
parallax measurement as well as investigate the functions which
now bear his name. Click HERE to
see a postage stamp issued in his honor showing Bessel and his
famous function.**

**GAMMA
FUNCTION: -The gamma function represents a
generalization of the factorial and is defined by the integral G(x)=Int[t^(x-1)*exp-t,{t,0,Inf}]. For
integer values of x one has that G(x+1)=x!
and
in general the generating function G(x+1)=x*G(x) holds. The function is encountered
in the series form for J(n,x) , whenever n is a
non-integer. One has G
(0.5)=sqrt(p)=1.77245 , G(1)=0!=1, and G(2)=1!=1.
A good table giving G(x) between x=1
and x=2 is all that is needed to find G(x)
at any other real value of x.**

**BESSEL
FUNCTIONS OF THE SECOND KIND [Y(n,x)]: -We show
here the first three Bessel functions of the second kind. Note
that they all go to minus infinity as x goes toward zero . This
is expected since x=0 is a regular singular point of the
equation. Although the Abel identity could be used to generate
this second linearly independent solution to Bessel's equation,
it is simpler for evaluation purposes to define it via the
indeterminate ratio Y(n,x)=lim n->integer{[cos(n*p )*J(n,x)-J(-n,x)]/sin(n* p)} due to Weber and Schlaefli and
sometimes referred to as the Neumann function. Note that if n is
a non-integer than the second linearly independent solution is
simply J(-n,x). The Hankel functions are the linear combinations
H1(n,x)=J(n,x)+i*Y(n,x) and H2(n,x)=J(n,x)-i*Y(n,x).**

**GENERATING
Y(0,X)
AS A LIMITING PROCEDURE: We know from the
Weber-Schlaefli indeterminate ratio that Y(0,x) can be generated
from the indeterminate ratio [cos(n* p)*J(n,x)-J(-n,x)]/sin(n*p) as n goes to an integer and that this
leads to rather complicated series expansions for the Bessel
Function of the Second Kind. (Go HERE
to see the rather lengthy derivation for Y(n,x) .) An
alternative approach for finding Y(n,x) is to take n very close
to the desired integer, in which case the above ratio is not
indeterminate and an evaluation involving the linearly
independent functions J(n,x) and J(-n,x) can be carried out. We
show you here the approximation to Y(0,x) obtained by letting
n=0.001. The resultant MAPLE output compares , as expected, very
well with the exact solution. Note that Y(0,x) becomes unbounded
at the singular point x=0 of the equation.**

**HYPERBOLIC(
OR MODIFIED) BESSEL FUNCTIONS: -The differential
equation x^2y"+xy'-(x^2+n^2)y=0 is known as the modified Bessel
equation and differs from the standard Bessel form only in that
a minus sign appears before the x square in the third term. The
simple substitution x-> - ix will map the two equations into
each other and thus the two lineraly independent solutions
I(n,x)=i^(-n)*J(n,i*x) and K(n,x) =( p/2)*[I(-n,x)-I(n,x)]]/sin(n* p) (known as the hyperbolic Bessel
functions of the first and second kind respectively) are
directly obtainable from the infinite series for J(n,x) . A
graph of the first few of these functions is shown in the graph.
Note that K(n,0) and I(n,inf) approach infinity and that one has
no zeros for finite x.**

**KELVIN FUNCTIONS: -These
functions are the solution of the generalized Bessel equation
x^2*y"+x*y'-(I*x^2+n^2)*y=0 and hence yield two solutions one of
which is y1=J(n,I^1.5*x). This function is complex with a real
and an imaginary part. Looking at the special case of n=0, one
finds that y1=Ber(x)+I*Bei(x), with Ber(x) and Bei(x) being the
convergent infinite series referred to as the Kelvin functions.
They arise in problems describing oscillatory viscous flow in
pipes and also in discussing the skin effect in high frequency
AC current flow through cylindrical conductors. Lord
Kelvin(family name- William Thompson) lived from 1824-1907
spending some 50 years as professor at the University of
Glasgow. Considered the greatest applied scientist and
mathematician of the Victorian era. He gave an estimation of the
age of the earth as 24 million years based on a heat transfer
calculation and his name is also linked with the laying of the
first TransAtlantic Telegraph cable in 1866. The earth's
age calculation was much too short(it's more like 3.7 billion
years) and not consistent with fossil evidence of life dating
back 200 million years and more. Click on Kelvinto
read
more about him.**

**SOLUTION
OF THE GENERALIZED BESSEL EQUATION: We have shown in class that the generalized
Bessel equation z^2y"+z(1+2m) y'+[m^2+((na)^2)z^(2n)-(nn)^2] y=0
has the general solution y=(1/z^m)[AJ( n,az^n)+BJ(-n,az^n],
where A and B are arbitrary constants. Let us use this result
to solve the problem y"+z^2y=0 subject to y(0)=1
and y (0)'=0 . Comparing this
equation with the generalized Bessel form one sees that
m=-1/2, n=+1/4 or -1/4, and a=1/2.
Thus the general solution is y=sqrt(z)*[AJ(1/4,z^2/2)+BJ(-1/4,z^2/2)]
and,
after application of the end conditions, one finds that A=0
and B= G (3/4)/sqrt(2). That is y=sqrt(z/2)* G(3/4)*J(-1/4,
z^2/2) . A MAPLE plot of this function is given in the
accompanying graph.**

**HYPERGEOMETRIC
EQUATION
OF GAUSS: -This equation reads
x*(1-x)*y"+[c-(a+b+1)*x]*y'-a*b*y=0, where a, b and c(or their
greek equivalents)are specified constants. The equation has
three regular singular points at x=0, 1 and infinity. The two
linearly independent solutions about x=0 are y1=F(a, b; c;
x)=1+a*b*x/c+a*(a+1)*b*(b+1)*x^2/[c*(c+1)*2!]+... and
y2=x^(1-c)*F(a-c+1, b-c+1; 2-c; x). Here F(a,b,c,x) is referred
to as the hypergeometric series and its radius of convergence is
one. Note that y2 will be identical with y1 when c=1 and hence
for this special condition the second solution needs to be
determined by use of the Abel identity. We show you here the
solutions y1=(1-x) and y2=1+(1-x)*ln[x/(x-1)] for the HGE when
a=c=1 and b=-1. Karl Friedrich Gauss(1777-1855) was a true
genius comparable to Archimedes and Newton and is often
referred to as the"prince of mathematicians". His interest were
wide ranging including, in addition to mathematics,
astronomy as director of the Goettingen observatory,
geomagnetism, geodesy and other areas of physics. Best known for
his proof of the fundamental theorem of algebra and things like
gaussian quadrature and the least squares method . If you click
HERE
you will see an image of him as it appeared on the German ten
mark bill. Also shown on this same bill is the gaussian
y=exp(-x^2).**

**LEGENDRE
POLYNOMIALS: - Graphs of
the first five Legendre polynomials Pn(x) , where P0=1, P1=x,
P2=(3*x^2-1)/2 , P3=(5*x^3-3*x)/2 and P4=(35*x^4-30*x^2+3)/8.
They represent the polynomial solutions to**

**[1- x^2]Pn"-2xPn+[n(n+1]Pn=0.
Adrien-Marie Legendre (1752-1833) was a member of the French
Academy of Sciences and well known for his mathematical
calulations including the gravitational ellipsoid problem
which led to the functions Pn(x). Click HERE to see an image of AML. He looks a lot like
the late Charles Laughton portraying Captain Blye in the movie
"Muntiny on the Bounty".**

**GENERATING FUNCTIONS FOR P(n,x): -The Legendre polynomials can be obtained via
several different generating functions. All of these will be
derived when we come to the section on integral
representations of the solutions of differential equations and
make use of the Euler kernel and the Schlaefli integral.
The inverse square root relation was historically the
first used to generate the P(n,x)s and follows directly from a
gravitational potential calculation for an asymmetric body. We
show you HERE the first nine Legendre polynomials
obtained via a simple one line program in MAPLE using
the generating formula P[n+1]={(x(2n+1)P[n]-nP[n-1]}/(n+1)
.**

**CHEBYSHEV
POLYNOMIALS: -These
represent one of the solutions of the Chebyshev equation
(1-x^2)*y"-x*y'+n^2*y=0 whenever n is an integer. Here are two
ways to generate their values from the differential equation.
The first approach is to introduce the independent variable
change t=(1/2)*(1-x) into the equation. This leads to a
hypergeometric equation with a solution y1(t)=F(n,-n;1/2;
t)=F[n,-n;1/2;0.5*(1-x)] which represents the Chebyshev
polynomials denoted in the literature by T(n,x). A second way
is to make the substitution x=cos(t) which leads to the
constant coefficient equation y(t)"+n^2*y(t)=0 which has the
obvious solution y(t)=cos(n*t). From this follows that
y1(x)=T(n,x)=cos[n*arccos(x)] and hence T(0,x)=1, T(1,x)=x,
T(2,x)=2*x^2-1, T(3,x)=4*x^3-3*x, and T(4,x)=8*x^4-8*x^2+1.
These polynomials are encountered in numerical analysis were
they are used to approximate functions within the range
-1<x<1. For example, the cubic y=x^3 equals
(1/4)*T(3,x)+(3/4)*T(1,x). The Chebyshev polynomials are
orthogonal to each other in abs[x]<1 for a weight
function of 1/sqrt(1-x^2). Furthermore, -1<T(n,x)<1 in
the same range of x as shown in the graph. Pafnuty Chebyshev
(1821-1894) was a professor at St.Petersburg University and is
known for his work on the prime number theorem, three bar
mechanical linkages, and his polynomials T(n,x). Click HERE to see the first few T(n,x) polynomials
generated by T[n+1]=2xT[n]-T[n-1].**

**LAGUERRE
POLYNOMIALS: -These are polynomials which arise
by solving the special form of the confluent hypergeometric
equation x*y"+(c-x)*y'-a*y=0 for c=1 and a= -n , where n is a
positive integer. The solution of this ODE( corresponding to the
Laguerre polynomials) is the Kummer series
y1=M(-n,1,x)=1-n*x+(-n)*(-n+1)*x^2/(1*2*2!)+. The standard
symbol for these polynomials is L(n,x) and one has L(0,x)=1,
L(1,x)=1-x, L(2,x)=(2-4*x+x^2)/2!, L(3,x)=(6-18*x+9*x^2-x^3)/3!
and L(4,x)=(24-96*x+72*x^2-16*x^3+x^4)/6!. The polynomials are
orthogonal in 0<x<infinity for a weight function of
exp(-x) and they play a central role in the solution of the
Schroedinger Equation for the hydrogen atom. Edmond
Laguerre(1834-1886) was a professor at the Ecole Polytechnique
in Paris who specialized in analysis and geometry. Today he is
mainly remembered for the L(n,x) polynomials through their
role in quantum mechanics. Click HERE
to see the first few computer generated Laguerre Polynomials
using the generation formula L[n+1]=(2n+1-x)L[n]-n^2L[n-1].**

**MATHIEU
FUNCTIONS: -As the last of our previously stated
12 classical ODE's we examine the Mathieu Equation
y"+[a-2q*cos(x)]y=0. This equation differs from the others in
that it has a periodic coefficient and thus it seems natural to
ask if there are not solutions to this equation which are also
periodic. Indeed as q goes to zero one has the periodic
solutions y1=cos[sqrt(a)*x] and y2=sin[sqrt(a)*x]. To find
periodic solutions for non-zero q one tries both even and odd
Fourier series expansions of period Pi or 2*Pi. This leads to
the four Mathieu functions ce(2*n+1,x), ce(2*n,x), se(2*n+1,x)
and se(2*n,x). In explicit form
ce(2*n+1,x)=A0*cos(x)+A1*cos(3*x)+A2*cos(5*x)+ where An
are coefficients determined by substituting into the
differential equation and using trignometric identities. Such
periodic solutions exist only for specified values of the
coefficients a and q as determined by evaluating the Hill's
determinant Det{Hill[a,q]}=0. In the graph we show you the
values of a and q corresponding to the ce(2*n+1,x) . Sorry for
the poor quality of the figure. This stems from the fact that
the Hill Determinant is available to me only via MATHEMATICA and
, as most of you know , MATHEMATICA is great for evaluation of
all sorts of functions and carrying out complicated evaluations
but is not very good when it comes to generating graphs.
You can also see a graph of the even and odd periodic fuctions
ce(2*n+1,x) and se(2*n+1,x) in the range of -2*Pi<x<2*Pi
by clicking HERE.
The value of q is 5 in both cases which requires a=1.85818754
and a=-5.79008060, respectively. Emile Mathieu (1835-1890) was a
French applied mathematician known for his work in group
theory, potential theory and for his 1868 paper on vibrating
elliptic membranes in which his equation and functions first
appear.**

**EXPONENTIAL
INTEGRAL VIA AN ASYMPTOTIC EXPANSION:-The
exponetial intergal is defined as
Ei(x)=Int[exp(t)/t,{t,-infinity,x}]. A simple integration by
parts allows one to expand it in the asymptotic form
Ei(x)=(1/x)*exp(x)*[1+1/x+2!/x^2+3!/x^3+O(1/x^4)] where the term
in the square bracket is a non-convergent asymptotic series
which will nevertheless give an answer for the value of the
integral close to its exact value when the number of terms
retained in the series equals the value of x. Thus, for example,
the Abramowitz and Stegun Math Handbook gives the
numerical value Ei(10)=2492.2289 while the
asymptotic form with ten terms yields 2492.49. The graph shown
plots the absolute value of the difference between Ei(10) and
the asymptotic series using n terms. It clearly confirms that
the best approximation occurs when n=x. You also can see a
comparison between Ei(x) and a seven term asymptotic
approximation by clicking HERE.**

**INTEGRAL
SOLUTION OF THE AIRY EQUATION USING A LAPLACE KERNEL: -
We have shown in class that second order ODE's with coefficients
linear in x can have an integral solution of the form
y(x)=const.Int[exp(xt)*v(t),{t,a,b}], where v satisfies the
first order ODE -d(Bv)/dt+C=0. Also one requires that the
bilinear concomitant P=vBK , where K=exp(x*t) is the Laplace
kernel, vanishes at the end point t=a and t=b. For Airy's
equation one finds that B=-1 and C=t^2. Thus the integral
representation for the Airy functions becomes Ai(x)(or
Bi(x))=Const*Int[exp(xt-t^3/3),{t,a,b}]. For large values
of complex t, the dominant term in P is
Real[exp(-t^3/3)]=exp-(1/3)*r^3*cos(3*theta) , where t
=r*exp(i*theta). So we see that P vanishes at large r in those
sectors of the t plane where cos(3*theta)>0. I show you these
sectors in the accompanying graph. It turns out that the Ai(x)
function is generated by going from a point 'a' at infinity in
the third quadrant to point 'b' at infinity in the second
quadrant of the t plane as indicated by the magenta colored
curve. The result of the contour integration connecting the two
points is found after a little manipulation to be
Ai(x)=(1/Pi)*Int[cos(x*u+u^3/3),{u,0,Inf}], where u=i*t. The
value of Const is established by looking at the integral in the
limit of vanishing x where the value should go to
Ai(0)=0.355028. If you click HERE
you will see the computer output for Ai(x) as obtained from its
integral representation using MAPLE. Note in obtaining this
result we have approximated the upper limit on the integral by
10 and thus one finds small wiggles for positive x which are not
present in the true Airy function.**

**INTEGRAL
REPRESENTATION OF J(n,x): -In the last few
lectures we have shown how one can obtain integral
representations to various functions which are governed by
second order ODEs with coefficients linear in x by use of a
Laplace kernel K=exp(x*t). Although the standard Bessel equation
is not in a form with linear coefficients, we can convert
it to one via the substitutions w=y*x^n and z=x^2. This leads to
the ODE z*w"+(1-n)*w'+(1/4)*y=0 which can be solved as
w(z)=Const*Int[exp(z*t-1/4*t),{t,alpha,beta}]. On converting
back to the x and y variables , this leads to the closed line
integral J(n,x)=[1/(2*Pi*x)] Int[exp(0.5*x(u-1/u)/u^(n+1)] about
the n+1 order pole at u=2*x*t=0. In turn, this line integral is
equivalent to the real integral representation J(n,x)=(1/Pi)*Int[cos[n*theta-x*sin(theta)]
,{theta,0,Pi}] as seen by
integrating about a unit radius circle in the u plane.
We show you here the result for J(0,x) generated
numerically from this last integral by use of MATHEMATICA.**

**STIRLING
APPROXIMATION FOR N FACTORIAL: -One of the
better applications of the Laplace method for determining the
asymptotic value of an integral with a large parameter deals
with evaluating n! for large n. Here one has the integral
n!=Int[t^n*exp(-t),{t,0,Inf}] whose integrand has one maximum at
t=n and thus n! has the approximate value (by the Laplace
method) of sqrt(2*pi*n)*n^n*exp(-n)
for n>>1. This is the famous Stirling approximation shown
in the graph. You will note that it already gives a very good
approximation to n! for values of n as low as 2. The Stirling
approximation finds wide application in the area of statistical
mechanics. James Stirling(1692-1770) was born and died in
Scotland, attended both Glasgow and Oxford University, and had
among his mathematician acquaintances Newton, Euler, N.Bernoulli
, de Moivre and McLaurin. Stirling was a mathematics teacher and
mining engineer but was unable to obtain a permanent university
position because of his support of the Jacobites. Stirling's
approximation first appears in his 1730 book Methodus
Differentialis , although de Moivre (a Huguenot expelled
from France and living in England ) may have already discovered
this formula a few years earlier. Go HEREto
see
a pdf file showing details of the derivation for the next term
in the asymptotic series for n!**

**AN
EXACT EVALUATION FOR N FACTORIAL: Although the
Stirling Approximation gives a pretty good approximation for the
factorial of large integers, one still very often requires exact
values for N! . Doing this by the brute force evaluation
of N!=product(n, n-1..N) can become rather time consuming.
Therefore one looks for some alternative approach to such an
evaluation. By clicking on the title of this section, you can
see one such attempt we developed for the evaluation of N! based
on an extention of the standard formula for (2N)! to (2^k*n)!.
We find, for example, the exact value of 64! using k=4 and n=4.
It is an integer with 89 digits and one your hand calculator
won't be able to display completely! Whether the present
approach is faster than the standard method for finding N! needs
to still be investigated and will depend on how much faster the
evaluation of 1.3.5.7...(n-1) =2^n*Gamma(n+1/2)/sqrt(Pi) is
compared to an evaluation of 1.2.3.4...n =n! when n gets large.**

**ASYMPTOTIC
FORM FOR THE BESSEL FUNCTIONS: -An application
of the method of stationary phase to Bessel's differential
equation leads to the asymptotic result
J(n,x)=sqrt[2/(pi*x)]*cos(x-pi*n/2-pi/4) for x>>1, with
Y(n,x) having same form except that cos is replaced by sin. We
plot here the result for J(0,x) and Y(0,x) and see that the
asymptotic forms are very good approximations to these Bessel
functions above x=2.**

**INTEGRATION
PATH
IN THE COMPLEX t PLANE FOR DETERMINING Ai(x) WHEN
x>>1: -We have shown in class that the
saddle point approximation to the integral
I=Int[Exp(x*h(t)),{t,a,b}] is given as
sqrt[2*Pi/(x*abs(h"))]*exp[x*h+i*(Pi-phi)/2], where phi is the
argument of the complex function h(t)" and evaluations are done
at the location of all saddle points t=tn along the path between
end points a and b. As is readily seen, the Laplace method and
the stationary phase technique are special cases of this more
general result. Applying the saddle point technique to the
integral form of the Airy function
Ai(x)=Const*Real{Int[exp[x*t-t^3/3],{t,-Infinity,+Infinity}]}
for large positive x one needs only to cross the saddle point
at t=-sqrt(x) to get from a=-Infinity to b=+Infinity. The
path which is taken is shown by the white parabolic curve in the
figure. The resultant approximation for the Airy function is
thus Ai(x)=(1/(2*sqrt(Pi))*x^(-0.25)*exp[-(2/3)x^1.5]
after adjustmwent of the constant. This result agrees closely
with the tabulated value of Ai(x) whenever x is greater than
about two. At x*=6.082202 the tables give Ai(x*)=8.101*10^(-6)
while the asymptotic approximation yields Ai(x*)=8.155*10^(-6).**

**ASYMPTOTIC
FORMS FOR THE AIRY FUNCTION VIA THE SADDLE POINT TECHNIQUE:
-In using the saddle point technique we found that for
x>0 one crosses the sadle at t=-sqrt(x) only in
determining the asymptotic form for Ai(x). When dealing with
negative values of x one needs to cross two saddles at
x=+i*sqrt[Abs(x)] and at x=-i*sqrt[Abs(x)] in going from t=a to
t=b. The results of such manipulations lead to the asymptotic
forms shown in the accompanying graph. Note that these
approximations are very good for x<-2 and x>2 but , as
expected, do not approximate Ai(x) well in -2<x<2.**

**COMPLETE
ELLIPTIC INTEGRALS OF THE FIRST AND SECOND KIND: -These
are intergrals encountered in determining the period of a simple
pendulum and in finding the circumference of an ellipse,
respectively. They are given by the infinite series
F(Pi/2,k)=K(k)=(Pi/2)*{1+[1/2]^2*k^2+[(1*3)/(2*4)]^2*k^4+
} and E(k)=(Pi/2)*{1-[1/2]^2*k^2-[(1*3)/(2*4)]^2*(k^4)/3-
}. Go HERE to see how
these infinite series expansions are obtained. The actual
numerical evaluation of complete elliptic integrals is most
efficiently accomplished by the alegbraic-geometric mean
(AGM)method of Gauss. Click HERE to see the
procedure.**

**JACOBI
ELLIPTIC FUNCTION Sn(u,k): -In determining the
angle versus time for the swing of a simple pendulum one
encounters the non-linear ODE (dx/du)^2=(1-x^2)*(1-k^2*x^2),
where u=sqrt(g/L)t , x=(1/k)*sin(theta/2)=Sn(u,k), t is the
time, and theta the swing angle of the pendulum. Solving this
equation with the intitial condition Sn(0,k)=0 leads to
the series solution
Sn(u,k)=x(u,k)=u-(1+k^2)*u^3/3!+(1+14*k^2+k^4)*u^5/5!+. We have
plotted a normalized version of this Jacobi Sine Elliptic
Function, namely Sn(u*K(m),m), for several different values
of m=k^2=[sin(thetamax/2)]^2, where thetamax is the
maximum swing angle of the pendulum. Here K(m) is the complete
Elliptic Integral of the First Kind. Carl Gustav Jacobi
(1804-1851) was a mathematics professor at Koenigsberg in East
Prussia and worked on differential equations and elliptic
functions during his rather short life.**

**JACOBI
ELLIPTIC FUNCTION Cn(u,k): -The
Jacobi
Cosine
Elliptic Functions are given by a solution of the ODE
(dx/du)^2=(1-x^2)*(1-k^2*x^2) subject to the initial condition
x(0)=1. The series solution reads
Cn(u,k)=x(u,k)=1-u^2/2!+(1+4*k^2)*u^4/4!+.. and a normalized
version , namely Cn(u*K(m),m), is plotted in the accompanying
gragh for several different values of m=k^2. Here K(m) is the
Complete Elliptic Integral of the First Kind.**

**ANALYTIC SOLUTIONS OF Y"=A*Y+B*Y^3 BY
ELLIPTIC FUNCTIONS:We have
shown in class that the Jacobi elliptic functions sn(u,k),
cn(u,k) and dn(u,k) have second derivatives which equal
certain non-linear expressions of the respective functions.
Thus, for example, sn"(u,k)=-(1+k^2)*sn(u,k)+2*k^2*sn(u,k)^3.
This suggests that one should be able to solve non-linear
differential equations of the form y"=A*y+B*y^3 by means of
elliptic functions. We demonstrate this possibility here by
examining the non-linear oscillator problem y"+y=a*y^3 subject
to y(0)=0, y'(0)=1. The boundary conditions suggest we compare
things with sn(u,k) and thus try a solution of the form
y=c*sn(b*x, k). Substituting this assumed form into the ODE
then requires, by comparison with the sn"(u,k) expression
given above and the boundary conditions, that A=-b^2*(1+k^2),
c*b=1, and B=2*(b*k/c)^2. If we now take A=-1 and B=a=0.494,
then k=0.9, b=0.7433 and c=1.3453. That is, the analytic
solution of the oscillator problem becomes
y=1.345*sn(0.7433*x,0.9). We have plotted this solution in the
accompanying graph and (as expected) things agree very well
with a numerical result obtained via a Runge-Kutta approach.**

**PHASE
PLANE PLOT FOR THE SIMPLE PENDULUM: - Here we show the phase plane trajectories
given by V^2/(2gl)=cos(theta)+Const. The separatrix between
oscillatory and non-oscillatory motion occurs at Const=1. The
governing non-linear differential equation for this problem is
x"+(g/l)sin(x)=0, where x equals the swing angle theta and "
represents the second time derivative.**

**PHASE PLANE SOLUTION FOR THE NON-LINEAR
PENDULUM WITH HIGH INITIAL ENERGY AND FINITE DAMPING:-One
of the more interesting examples of a solution to a non-linear
differential equation is that for a simple pendulum with high
initial rotational speed and damping. For this case one has the
equation x"+ ax'+sin(x)=0
subject to x(0)=a and x'(0)=b. We have carried out a
Runge-Kutta numerical evaluation of this autonomous equation for
the case of a=0.3 and the ICs of a=0
and b=4. The two-line mathematical program using MAPLE and
the resultant phase trajectory (after enhancement via
paintbrush) is shown in the accompanying figure. The results
clearly indicate that the initial kinetic energy of the pendulum
is larger than that required to take it over the top, but that
the damping dissipates this energy in time and eventually the
pendulum goes into a damped periodic motion and finally comes to
rest at the critical point X=2p, X'=0
as the time approaches infinity. This numerical procedure
using MAPLE is simpler and faster than that required by either
MATLAB or MATHEMATICA, however, I am not very happy with the
convoluted way MAPLE expresses higher derivatives and writes out
ODEs.**

**PHASE PLANE
SOLUTION IN THE VICINITY OF CRITICAL POINTS: -We have shown that an autonomous non-linear
equation of the form F(x",x',x)=0 can be re-written as the
first order phase plane equation dy/dx=P(x,y)/Q(x,y), where
y=dx/dt.**

**Now there generally are one or more
points within the phase-plane where dy/dx becomes
indeterminate. These are the critical points of the equation
and one can linearize the functions P and Q near these points
to get a good idea of the type of phase trajectory expected in
the neigborhood by looking at the solution of the linear
fractional equation dy/dx=(a*x+b*y)/(c*x+d*y) whenever the
critical point lies at x=y=0.(A obvious shifts in coordinates
is required when the critical point is located at other than
the origin). Here a=dP/dx, b=dP/dy ,c=dQ/dx, and d=dQ/dy are
constants evaluated at the critical point. This equation can
be solved in closed form and the type of phase trajectory can
be determined by looking at the roots of (lambda)^2-(b+c)*(lambda)+(b*c-a*d)=0 . Real and opposite signs in lambda
correspond to a saddle point, while complex conjugates give
spirals , real lambdas of the same sign corrspond to a node,
and pure imaginary roots of opposite sign yield a center. To
support these observations we compare the exact phase
trajectories of x"+x+x^2=0 with the results predicted by
a linearization of dy/dx=-x*(1+x)/y . This last equation has
critical points at (0,0) and (-1,0) and substituting into the
lambda formula clearly shows the first is a center and the
second a saddle point. These results are seen to agree well
with the exact solution in the immediate neighborhood of the
two points. The solutions are Liapunov stable in time if the
real parts of the lambdas are negative, and unstable if
one or the other of the real parts is positive. The trajectory
about a center is neutrally stable.**

**PHASE PLANE
SOLUTION OF X"=X*(1-X^2): -This is another
interesting autonomous equation which is equivalent to the first
order ODE dy/dx=x*(1-x^2)/y and has the analytic solution
(1/2)*x^2-(1/4)*x^4-(1/2)*y^2=Const. A plot of this solution is
shown on the accompanying graph for several different values of
the constant. Note the bowtie configuration of the separatrix
and the location of the critical points at (-1,0) and (1,0)
representing centers and the one at (0,0) representing a saddle
point.**

**By just changing the sign on the right hand side of this phase
plane equation, one obtains the completely different trajectory
shown HERE which exhibits two saddle points and one
center.**

**A
NON-LINEAR EQUATION WITH NINE CRITICAL POINTS: -In
discussing
autonomous
differential equations we have studied the behaviour of their
phase plane solutions in the vicinity of critical points. One
can also work things backwards by first constructing a first
order phase plane equation with specified locations of the
critical points and then obtaining the corresponding second
order non-linear equation. One such manipulation starts with
dy/dx=[x*(1-x^2)]/[y*(1-y^2)] which clearly has nine critical
points corresponding to dy/dx=0/0 and which has the analytic
solution y^2*(2-y^2)-x^2*(2-x^2)=Const. A plot of this solution
is shown in the accompanying graph and the corresponding second
order ODE is x"=x(1-x^2)/(1-x'^2). Note that there are 4 centers
and 5 saddle points. The solutions are periodic whenver the ICs
lie within the circle x^2+y^2=2 but are unbounded for any
starting point outside this circle. The implicit solution t=t(x)
can readily be found by an exact integration of the phase plane
trajectories and is t=A+Int[1/sqrt[1+sqrt(B-x^2(1-x^2))]], with
A and B being constants determined by the ICs x(0) and x'(0). A
numerical integration of this last integral for the starting
conditions x(0)=0.5 and y'(0)=0 shows the solution has a period
of tau=1.105384..**

**VOLTERRA
PREDATOR-PREY PROBLEM: -A
problem which is well treated by phase plane techniques is the
predator-prey problem of Volterra in which the interaction
between prey y and predator x are modelled according to
the simultaneous equations dy/dt=(a-b*x)*y and
dx/dt=(-c+d*y)*x, where a, b, c, and d are constants. Dividing
the first of these by the second yields a first order phase
plane ODE which has the analytic solution
b*x+d*y-log[(y^c)*(x^a)]=Const. One can either plot this
solution directly to get a phase plane plot or attack the
original set of simultaneous first order equations directly
using the Runge-Kutta numerical method. The plot we show here
was gotten using the numerical approach via MAPLE. From the
plot one sees that there is a center at x=a/b and y=c/d and
that the population sizes of both x and y are periodic
functions of time t. Vito Volterra (1860-1940) was a
mathematics professor at the Universities of Pisa, later Turin
and then Rome specializing in integral equations and
population growth problems for interacting species. He lost
his teaching job in 1931 for refusing to sign an oath to the
Mussolini regime.**

**LIMIT
CYCLES: Certain nonlinear second order
equations exhibit the interesting property that their phase
plane trajectories merge to a single closed curve away from any
critical points of the problem as time becomes large. Such
trajectories are referred to as limit cycles and , according to
the Poincare-Bendixson Theorem, are located somewhere in the
annular region between two concentric circles in the phase plane
if all the phase trajectories enter this annular region from the
outside with increasing time. Note that most non-linear
autonomous equations do not have limit cycles. We show you here
a numerical solution for one special nonlinear equation which
happens to have x^2+y^2=1 as its limit cycyle.**

**TRANSIT
TIME ABOUT A LIMIT CYCLE:
It is often of interest to find the time it takes the solution
within the phase plane to make a complete transit around the
limit cycle as this will give the period of the solution x(t)
at large time. Such a determination is generally not possible
by analytic means, however, can be gotten by numerical
approaches. We demonstrate the procedure here by looking at
the equation x"+(x^6-1)x'+x=0. This equation clearly has a
limit cycle because the damping term x^6-1 is negative for
small x but positive for large x. The phase plane equations in
this case are dx/dt=y, and dy/dt=(1-x^6)y-x. We solve these
numerically by starting at t=0 at an arbitraryly chosen
initial point say x0=y0=0.2 . Running things for large
enough t will then locate the limit cycle for the problem.
Next we repeat the numerical solution by choosing t=0 at a
starting point x1,y1 somewhere on the limit cycle. Then note
the time tmax it takes for your solution to make one
loop around the limit cycle. This time will equal the
limit cycle period which here turns out to be tmax=7.34 .**

**VAN DER POL EQUATION AND THE LIMIT CYCLE: -This
is one of the first non-linear differential equations studied
historically which exhibits a limit cycle. Here is a numerically
determined phase plane plot for the VDP at e=0.1. Note,
regardless of the initial conditions, the solution eventually
ends up on the limit cycle, which for this small value of e, is
nearly a circle. Van der Pol used these results to explain
parametric oscillations observed in a triode oscillator. Note
that the Poincare-Bendixson theorem indicates that there exists
a limit cycle for this case lying within the annulus
0.5<sqrt[x^2+y^2]<2.5. A nice applet showing the phase
plane plot for the van der Pol equation for differernt epsilons
and ICscan be found by going to- http://www.apmaths.uwo.ca/~bfraser/version1/vanderpol.html
I thank one of our
students, Inka Hublitz, for calling my attention to this URL.
Note that the applet does not give a one to one projection, so
that what should be a circular limit cycle for very small
epsilon is projected as an ellipse.**

**RAYLEIGH
EQUATION: -Another
differential equation exhibiting a limit cycle is that of Lord
Rayleigh, namely, X"- m
(1-X'^2)X'+X=0. This equation has similarities with the van
der Pol equation and was formulated by Rayleigh in his Theory
of Sound book some thirty years before van der Pol. It
describes nonlinear velocity dependent damping in a
simple mass-spring system and is also encountered in the
mechanical problem of stick-slip motion. The equation has a
critical point at X=X'=0 and the trajectory within its
neighborhood is either an unstable spiral or an unstable node
depending on the magnitude of m.
We show you here its solution when m
=1 for two different initial conditions . Note that after a
short time, the solution ends up on the oval shaped limit
cycle regardless of where things started at t=0. This will
always be the case when a limit cycle exits and guarantees a
periodic behaviour of the solution X(t) once the limit cycle
is reached.**

**BOGOLIUBOFF-KRYLOFF
APPROXIMATION:
-This is an analytical approximation method for the
solution of non-linear equations of the type x"+e*g(x',x)+x=0,
where e(or more commonly epsilon) is a small parameter and
g(x',x) is a nonlinear function. Noting that for e=0 this
equation has solution x=A*sin(t+p) , where p is the constant
phase and A the constant amplitude, B&K assumed that for
small e the solution should have the almost periodic form
x=A(t)*sin(t+p(t)). Differentiating this expression with respect
to t then yields x'=a'*cos(t+p) provided one sets the remaining
portion of the differentiation, ,namely,
A'*sin(t+p)+A*p'*cos(t+p)=0. Differentiating once more and
substituting into the original ODE, one can solve for the
time derivatives of A and p. Averaging these over one period
leads to the BK approximation A'=-e/(2*Pi)*Int[g*cos(psi),{psi,
0, 2*Pi}] and p'=e/(A*2*Pi)*Int[g*sin(psi) ,{psi, 0, 2*Pi}),
where psi=t+p(t). Applying this approximation to the van der Pol
equation for e =0.1 yields
x(t)=A*exp[0.5*e*t)*sin(t+p(0))/sqrt[1+0.25*A^2*exp(e*t)]. In
the phase plane this yields the graph shown with a limit cycle
corresponding to an amplitude of A=2. The result is seen to be
very close to the numerical solution shown earlier.**

**RESPONSE
OF
A HARD-SPRING OCILLATOR TO A PERIODIC DRIVING FORCE: -A
method
to
solve non-linear equations of the Duffing type for a periodic
driving force F*cos(w*t) is to try a harmonic expansion of the
type x=A1*cos(w*t)+A2*cos(3*w*t)+A3*cos(5*w*t)+...For the
ODE x"+k*x'+a*x+b*x^3=B*cos(w*t), this leads, after some
manipulations involving use of various trignometric identities,
to the result (A1/B)^2=1/[(a-w^2+0.75*b*A1^2)^2+(k*w)^2]. This
response function is plotted in the accompanying graph. Note
that the resonance point is moved to the right of that
corresponding to the linear oscillator case at b=0. Its form
also explains the sudden change in amplitude response observed
experimentally as the driving frequency is increased
progressively from a points either above or below the linerar
resonance value at w/sqrt(a)=1.**

**SOLUTION OF
DUFFING'S EQUATION: -Here we show a numerical
solution to the non-linear Duffing equation when the system is
subjected to a sinusoidal input. The phase plane plot for the
special initial conditions chosen demonstrate an interesting
double attractor pattern and the displacement x(t) versus time
t exhibits a constant amplitude periodic response despite
the fact that a damping term is present in the equation.**

**STICK-SLIP MOTION: -Another area
where one encounters a phenomenon involving a limit cycle is the
stick-slip motion of a mass sitting an a constant velocity(v)
conveyor belt and connected via a linear spring(k) to a fixed
point. By expanding the resultant friction force F(dx/dt-v) in a
Taylor series about v and noting that friction force always
opposes the mass displacement x(t), one finds that
x"-ax'+b*(x')^3+x=0 approximately, where x now equals x+f(v)/k.
We show the phase plane solution for this equation in the
accompanying graph. Note near the origin the trajectories are
unstable spirals while for large values of y=x' one has an
inward moving trajectory. Hence from the Poincare-Bendixson
theorem, one can conclude that there is a limit cycle as the
graph shows.**

**MATLAB
ODE23 PROGRAM FOR SOLVING NTH ORDER NON-LINEAR ODES: -We
have
now
reached a point in the course where we are considering
non-linear and non-autonomous ODES. These no longer lend
themselves to most of the analytical solution techniques which
have been discussed. Although one can continue to try to solve
them via a Taylor series expansion, it is not at all clear what
the range of converence of the resultant series will be since
one does not know beforehand at what point the solution will
become unbounded. Also for non-linear equations, the points
where the solutions exibit a singularity is very much a function
of the imposed initial conditions . Under such
circumstances it is often better to just attack the
equation directly by a numerical procedure. One very convenient
way for doing so is to use the MATLAB canned program ODE23. By
studying the scan I have made for solving the sample equation
x"+t*sin(x)=0 subjected to x(0)=1 and x'(0)=0, you will see how
quickly such a numerical solution can be obtained. The ODE23
program uses the famous Runge-Kutta numerical stepping method
first introduced in 1901 by Carl Runge(1856-1927) and Martin
Kutta(1867-1944). Runge was a professor at Goettingen also
working on atomic spectra, while Kutta was a professor at
Stuttgart also known for the "Kutta Condition" that flow
leaves the trailing edge of an airfoil at finite speed.**

**EMDEN'S EQUATION: -About 1906
R.Emden was studying the radial density variation of stars based
on Newton's gravitational potential theory and the classical gas
laws for an ideal gas. In his studies he came up with the
equation x"+(2/t)*x'+x^n=0 which is now referred to as Emden's
equation. Here x is a measure of the gravitational potential and
t a non-dimensional radial direction. The constant
n=1/(gamma-1), where gamma is the ratio of specific heats of the
gas constituting the star. For a monotonic gas one has
n=1/(1.6666-1)=1.5. We show via MATLAB ODE23 the x(t) dependence
for n=1, 1.5, and 2. Note that for n=1 the equation becomes
linear and has a solution which just equals j(0,t)=sin(t)/t.**

**BLASIUS
EQUATION: -This is
the third order non-linear equation
x*(d^2x/dt^2)+2*(d^3xdt^3)=0 arising in connection with
viscous flow past a flat plate. First solved by H.Blasius
in
his
dissertation . Blasius was a student of L. Prandtl of boundary
layer fame. The story told to me by K. Pohlhausen is that
Blasius was so disappointed with his lack of progress toward his
PhD at Goettingen that he was ready to leave the University. It
was at that point that Prandtl grabbed him and told him to get
busy on solving a certain new higher order differential equation
which he (Prandtl) had discovered but couln't solve. Blasius
produced a solution using an analytic procedure where one
couples a Taylor series expansion for small t with an asymptotic
solution for large t . I show you here the MATLAB solution
again produced by the " garden hose" technique of guessing
a value of d^2x(0)/dt^2 until the solution meets the required
derivative condition at infinity. The value of this
constant is found to be 0.332 and can be used to calculate
the viscous drag a flat plate experiences under laminar flow
conditions.**

**THOMAS-FERMI
EQUATION: -Another interesting non-linear ODE of
the non-autonomous type is that of Thomas and Fermi. It reads
y"=y^1.5/sqrt(x) and arises in describing (by means of Poisson's
equation) the spherically-symmetric charge distribution about a
many electron atom. We give here a Runge-Kutta numerical
solution to the problem based on the "garden hose" approach. One
is using the physically meaningful boundary conditions y(0)=1
and y(inf)=0. It is of particular historical interest that this
was the first non-linear differential equation attacked by a
forerunner of the electronic computer , namely the 1931 MIT
differential analyzer of Bush and Caldwell . We point out that
the first true electronic computer was built a few years
later(about 1939) by John Vincent Atanasoff
, who had been an undergraduate here at the University of
Florida (BS in electrical engineering 1925).**

**NUMERICAL
SOLUTION
OF A NON-LINEAR ODE USING MAPLE: An even simpler to use canned program for
solving an nth order ODE numerically is provided by MAPLE. I
show you here the couple of lines required to solve and then
plot the phase plane solution and the x(t) solution of the
non-linear and non-autonomous ODE x"=-t^2*sin(x) subjected to
the ICs x(0)=1, x'(0)=0. Note the nearly periodic character of
the solution as t becomes large. The canned technique
employeed here is Runge-Kutta. To treat boundary value
problems using this approach one adjusts the left derivative
condition by trial and error until the solution falls through
the right boundary value point(ie-garden hose method).**

**SOLUTION
OF A NON-LINEAR BOUNDARY VALUE PROBLEM BY THE GALERKIN
METHOD-To complete our discussion on ODEs, we
consider solving a non-linear ODE by the Galerkin method. This
analytical approach works especially well for boundary value
problems and , unlike a numerical approach based on the Runge
Kutta technique, does not require knowledge of both x(0) and
x'(0). The basic idea behind the Galerkin approach is to
expand the equation solution as**

**x(t)=Sum[c _{n}f_{n}(t),
n=1..N], where f_{n}(t) are trial functions each of
which satisfy the boundary conditions of the problem. If one
now takes the governing equation x"(t)=F[x'(t),x(t),t]
multiplies it by f_{m}(t) and integrates over the
range 0<t<T, this will lead to a set of non-linear
algebraic equations for the c_{n}s which when solved
give an analytic approximation for x(t). We demonstrate this
approach for x"=x^{2 }with x(0)=x(1)=0. Choosing just
a single trial function f_{1}(t)=c_{1}sin(pt) , one finds c_{1}=_{ }-(3/8)p^{3}. That is, x(t)=-11.63sin(pt) . As seen from the graph this
result is remarkably accurate and could be further improved by
taking more terms in the Galerkin series. The technique also
works very well for the determination of eigenvalues in linear
differential equations.**

**EGM 6322-PARTIAL
DIFFERENTIAL EQUATIONS AND BOUNDARY VALUE PROBLEMS**

*Partial
differential
equations of first and second order. Hyperbolic, parabolic,
and elliptic equations including the wave, diffusion, and
Laplace equations. Integral and similarity transforms.
Boundary value problems* *of the
Dirichlet and Neumann type. Green's functions, conformal
mapping techniques, and spherical harmonics.*
*Poisson,
Helmholtz,
and Schroedinger equations.*

**SOLUTION
SURFACE FOR A FIRST ORDER PDE: -Exact analytic solutions exist for many
linear first order partial differential equations. Here is the
solution surface for one such equation, namely, dz/dx+dz/dy=xy
when z(x,0)=x^3. The solution reads
z=-(1/6)*(x^3)+(1/2)*(x^2)*y-(7/6)*(y-x)^3 and has the shape
of a flying carpet as shown in the figure.**

**SOLUTION OF A
FIRST ORDER PDE CONTAINING A SPACE CURVE: - We have given a general solution for the
linear first order PDE's of the form
f(x,y)dz/dx+g(x,y)dz/dy=0 in class and showed how this
solution can be adjusted to contain a single space curve. The
graph gives one such solution surface for ydz/dx-xdz/dy=0
containing the corkscrew curve x=t*sin(t), y=t*cos(t), and
z=t.**

**CRYSTAL GROWTH EQUATION: -The
growth of crystals in a crystalizer can be well described by
the first order PDE y _{t}+Ry_{x} +y/T=0 ,
where y is the crystal number density, x the crystal diameter
and T the drainage time for the fluid throughput in the
crystalizer. We show you here the solution
y=exp[-0.1t-(x-2*t)^2/20] , obtained by the characteristic
method, for three different times t , assuming an initial
gaussian distribution of y(x,0)=exp(-x^2/20). This type
of equation also is applicable to the problem of urinary
concretion growth (ie-kidney stone formation).**

**CHARACTERISTICS
OF
SECOND ORDER PDE'S:We have
shown that the 2nd order linerar PDE a(x,y)z _{xx}+2b(x,y)z_{xy}+c(x,y)z_{yy}=0,
with
variable coefficients a, b and c , can be cast into a
canonical form containing only one second derivative
when using the equation's characteristics eta and zeta as the
independent variables. To determine these characteristics one
needs to solve the first order ODE dy/dx=[bï¿½
sqrt(b^2-ac)]/a. For a derivation of the general
canonical form click HERE. Consider
next the special case of the Tricomi equation z_{xx}+xz_{yy}=0
where
a=1, b=0, and c=x. This equation arises in the study of
transonic flows. A simple integration yields the
characteristics h=y+(-x)^1.5
and z=y-(-x)^1.5
as shown in the graph. Note that this equation is hyperbolic
with real characteristics whenever x<0 and elliptic
with non-real characteristics for x>0 as determined by the
sign of the discriminant (b^2-ac)=-x.**

**CHARACTERISTICS OF THE PRANDTL-GLAUERT
EQUATION:An interesting
second order constant coefficient PDE is the Prandtl-Glauert
equation [M ^{2}-1]j_{xx} -j_{yy}=0, where M=U/c is the Mach number and j(x,y) the
velocity potential for a linearized version of the
steady-state 2D Euler equation combined with the divergence
and irrationality conditions for an inviscid compressible
flow. The slope of the characteristics are readily seen
to be tan(q)=dy/dx= ï¿½1/sqrt[M^{2}-1].
That is, the characteristic lines make an angle q=ï¿½ arcsin(1/M) with respect to the x axis. These directions
are equivalent to the famous Mach lines and correspond to the
angle a shock wave makes about a body moving at supersonic
speeds. Click on the above title to see a photo of such an
oblique shock as I was able to obtain via schlieren
observation of flow past a sharpened nail held in the test
section of a small supersonic wind tunnel which I constructed
for my high school science project many years ago. The angle
in the photo shows that air is moving past the projectile at
M=1.76. Such speeds are comparable with those of a high
powered rifle bullet or the Concorde supersonic airliner. Note
that the above PDE remains hyperbolic and hence has real
characteristics only as long as M>1.**

**D'ALEMBERTS
SOLUTION TO THE 1D WAVE EQUATION: -In his famous
paper of 1747, 30 year old Jean D'Alembert showed that the
solution for a vibrating string of infinite length is
z(x,t)=(1/2)*[f(x-ct)+f(x+ct)]+(1/2c)*Int[g(t),{t,x-ct,x+ct}],
where z(x,0)=f(x) and dz(x,0)/dt=g(x). Here we demonstrate this
result for the special case of a gaussian initial displacement
z(x,0)=exp-x^2 and no initial velocity g(x)=0. Jean Le Rond
D'Alembert(1717-1783) was the illegitimate son of Mme de Tencin
and an artillery officer and was abandoned as an infant on the
doorsteps of the Le Rond church in Paris. A very brilliant
but quarrelsome individual (he had disputes with Bernoulli,
Clairaut, and Euler among others), D'Alembert worked most of his
life for the Paris and French Academy of Sciences and was also a
member of the Royal Society of England and of the Prussian
Academy of Sciences . Besides the vibrating string problem, he
is best known for the D'Alembert Principle in mechanics.**

**FOURIER
TRANSFORM FOR THE RECTANGULAR PULSE: -The
Fourier transform of a function f(t) is defined as
F(w)=Int[f(t)*exp(-iwt),{t,-Inf,+Inf}]. Its inverse is
f(t)=(1/2Pi)*Int[F(w)*exp(iwt),{k,- Inf,+Inf}]. We show here
f(t) and its transform F(w) for the rectangular pulse f(t)=+1
for -1<t<1 and f(t)=0 for all other t. Note that in
physics one uses a slightly different and more symmetric
definition of the Fourier transform in which the constant in
front of both integrals is 1/sqrt(2Pi) and the sign on the
exponentials is changed. It makes no difference which form is
used as long as one remains consistent. Joseph
Fourier(1768-1830) taught at the Ecole Polytechnique,
accompanied Napoleon on his Eqyptian campaign as scientific
advisor, and was appointed prefect of Grenoble. Best known for
his book on the theory of heat conduction in which he developed
his famous Fourier series expansion.**

**FOURIER INTEGRAL SOLUTION OF THE 1D WAVE
EQUATION: -An alternative way to solve the 1D
wave equation in the infinte x range is to use a Fourier
Transform of the form
F[f(x)]=Int[f(x)*exp(-i*k*x),x=-infinity..+infinity]. This leads
to an ODE which can be solved and then inverted to recover the
explicit solution U(x,t). We show here such an approach for a
string of infinite length subjected to an initial rectangular
pulse displacement and no initial velocity.**

**FOURIER SERIES FOR THE TRIANGLE FUNCTION: -This
shows
a 31 term Fourier Series approximation to the even periodic
function F(x)=Abs[1-x] in -1<x<1. There is no Gibbs
overshoot in the series representation of F(x) at integer values
of x since there is no jump in the value of the function at its
slope discontinuities.**

**GIBBS
PHENOMENON: -This represents the overshoot
observed in a Fourier series expansion of a function at a
discontinuity of that function. Here we show this phenomenon for
the rectangular pulse F(x)=+1 for 0<x<Pi and F(x)=-1 for
-Pi<x<0 and magnify the behaviour of a 100 term Fourier
series approximation near the discontinuity at x=0. This
overshoot was first observed by Albert Michelson(of
Michelson-Morley experiment and speed of light fame) in
his Fourier series generator and first explained by vector
field expert Willard Gibbs of Yale . The overshoot does not
dissapear even as the number of terms approaches infinity,
however the area under the overshoot does approach zero.
According to Dirichlet, a Fourier series converges to a point
representing the mean value of the function F(x) on the two
sides of a discontinuity.**

**SEPARATION OF VARIABLES SOLUTION FOR A
VIBRATING STRING:We have
shown in class that the 1D wave equation in a finite x region
is simplest to solve via the separation of variables approch
leading to a Fourier series expansion. In general, for a
solution in 0<x<L with boundary conditions
U(0,t)=U(L,t)=0, one finds that U(x,t)=Sum[sin(n*Pi*x/L)*(a _{
n}*cos(n*Pi*c*t/L)+b_{ n}*sin(n*Pi*c*t/L)),
{n=1, infinity}], where a_{n}=(2/L)*Int[U(x,0)*sin(n*Pi*x/L),{x,0,L}]
and b_{n} =(2/L)*Int[U(x,0)*cos(n*Pi*x/L),{x,0,L}]. We
show here the solution (obtained via MAPLE) for a vibrating
string with an initial triangular displacement U(x,0)=x for
0<x<0.5 and U(x,0)= (1-x) for 0.5<x<1 . The
propagation speed has been set to c=1 , we used the first 60
terms in the Fourier expansion, and set the time
at 1/8 th period intervals.**

**EIGENMODE FOR
A RECTANGULAR MEMBRANE: -The
wave equation for a vibrating rectangular membrane
clamped at its edges is governed by U _{tt}=c^2[U_{
xx}+U_{yy}] subjected to the four boundary
conditions U(0,y,t)=U(a,y,t)=U(x,0,t)-U(x,b,t)=0 plus two
initial conditions U(x,y,0)=phi(x,y) and U_{t}(x,y,0)=psi(x,y).
A
simple separation of variables approach gives the double
Fourier series solution taken over the integers n amd m. The
spatial portion of this solution is composed of the eigenmodes**

**VIBRATING
CIRCULAR MEMBRANE: -The wave equation yields
standing wave solutions on a circular membrane clamped at its
edge which are expressible as the product of radially dependent
Bessel functions and an angularly dependent cosine term. Each of
the fundamental modes has associated with a unique oscillation
frequency. A contour plot of the wave pattern associated with
an n=3 and m=2 mode for a unit radius membrane is shown in
the figure. Also an oblique view of an axisymetric mode is found
by clicking HERE.**

**IMPLODING
GAUSSIAN: The 3D wave
equation for a spherically symmetric wave reads U _{tt}=c^2*[U_{rr}+(2/r)*U_{r}
]. The substitution U=V(r,t)/r leads to an equivalent 1D wave
equation V_{ tt} =c^2*V_{ xx}, which has the
well known D'Alembert solution. A special case of this
solution is a spherically symmetric imploding wave given by
U(r,t)=(1/r)*F(r+c*t) corresponding to an initial displacement
of U(r,0)=F(r) and zero initial velocity U_{t}(r,0)=0.
We show here such an imploding wave when the initial
displacement is the gaussian F(r)=exp[-100*(r-1)^2].
Note the rapid increase in wave amplitude as the wave travels
toward the origin. This concept finds application in
nuclear bomb triggers and is also connected to the problem of
laser induced deuterium pellet ignition in nuclear fusion.**

**SOLUTION
OF A NON-HOMOGENEOUS WAVE EQUATION: A variation on the solution of the standard
wave equation arises when the equation contains an extra term
F(x,t) so that things read V _{tt}=c^{2}V_{xx}+F(x,t)
subject
to say V(0,t)=V(a,t)=0 and V(x,0)=V_{t}(x,0)=0.
In this case one still has V(x,t)=Sum[C(t)*sin(np/a), n=1..infinity], but since one can
expand F(x,t) as Sum[D(t)*sin(np/a),
n=1..infinity], one needs to satisfy the ODE**

**TEMPERATURE
IN A 1D BAR HAVING AN INITIAL DELTA FUNCTION INPUT: -Using the general solution of the 1D
temperature equation in a bar extending from minus to plus
infinity, we find the following graph for an initial delta
function input at t=0 . The plots are for three different
times(thermal diffussivity has been set to unity) and since
there is no heat loss in the process the total energy and
hence the area underneath these curves stays constant despite
of the diffusion which is occuring. The solution T(x,t) can be
considered as a continous representation for the delta
function as t is allowed to approach zero.**

**SOLUTION OF
THE 1D DIFFUSION EQUATION: -We have shown via a
Fourier transform that the one dimensional diffusion equation in
the infinite spatial range -Inf<x<+Inf has the integral
solution
C(x,t)=[1/(2*sqrt(D*t*Pi)]*Int[f(z)*exp-[(x-z)^2/(4*D*t)],{z,-Inf,+Inf}].
Here
the initial condition is C(x,0)=f(x) , z is a dummy variable,
and D is the constant diffusion coefficient. Making use of
symmetry, one can use this result to solve the problem of
diffusion of a substance in the half-infinite range
0<x<Inf for the initial distribution C(x,0)=1 in
0<x<1 and the non-penetration condition dC(0,t)/dx=0 and
vanishing condition C(Inf,t)=0. The result is shown in the
accompanying graph. Note that about one-half of the intitial
concentration has dispersed from the initial non-zero region at
time t=a^2/D, where a=1 is the initial concentration width.**

**ERROR FUNCTION
[ERF(X)] AND ITS COMPLEMENT [ERFC(X)]: -This graph shows the value of the
integral [2/sqrt(pi)]*Int[exp-x^2,{t,0,x}] encountered in
studying the heat flow into a semi-infinite bar of zero initial
temperature subjected to a constant temperature at its left end.
The error function is also encoutered in certain diffusion
problems such as in determining the spatially dependent doping
distribution in semi-conductors.**

**APPLICATION OF DUHAMEL'S PRINCIPLE TO SOLVE
THE 1D HEAT CONDUTION EQUATION :Jean Marie Constant DuHamel(1797-1872) , who
was a professor at the Ecole Polytechnique in Paris,
introduced a technique which allows one to express the
solution of the 1D heat conduction equation with
time-dependent end conditions in terms of the much simpler
known solution when the end condition is constant. The
technique, now referred to as the DuHamel principle, is
based on the use of Laplace transforms. Lets demonstrate for
the problem T _{t}=T_{xx }in
0<x<infinity if T(x,0)=0, T(0,t)=sin(t) and
T(infinity,t)=0. Taking the Laplace transform with respect to
time and applying the initial and the two boundary conditions
, one finds T_{bar}=[s/(s^2+1)]*[(1/s)*exp-(sqrt(s)*x)].
But
we recognize that the terms in the two square brackets are the
transforms of just cos(t) and erfc(0.5*x/sqrt(t)),
respectively. Hence, by the convolution theorem, one has that
T(x,t)=Int[cos(t-v)*erfc(0.5*x/sqrt(v),{v,0,t}]. We show you
here a plot of this function at t=Pi by numerically evaluating
the integral with aid of MAPLE. It takes some 400 seconds of
calculation time on my machine to obtain this single curve.**

**THERMAL WAVES:Although the heat
conduction equation does not yield a wave solution in the
standard form, it does allow the existence of a highly damped
and dispersive temperature signal into a conductor when the
surface temperature is varied periodically. We demonstrate this
point by looking at the solution of T _{ t}=aT_{xx} for 0<x<infinity
when T(0,t)=sin(w*t) and
T(infinity,t)=0. At larger time( where initial temperature
conditions are no longer important) one can assume that
T(x,t)=Imaginary[exp(i* w*t)*F(x)] .
This leads to the solution T(x,t)=exp(-sqrt[
w/(2*a )]*x)*sin[w*t-sqrt( w/2a )*x] which we have plotted in the
accompanying graph for the case of copper where a=1.13cm^2/sec at four different times.
The average penetration of the highly damped thermal wave is
seen to be approximartely equal to x(in cm)=sqrt[a/t] as is to
be expected from pure dimensional arguments. Such time dependent
temperature variations can be used to make thermal diffusivity
measurements in thin films using laser generated heat pulses.**

**SOLUTION OF THE HEAT CONDUCTION EQUATION FOR
A FINITE LENGTH BAR: -The
solution to the heat conduction equation within a finite
length bar follows readily by a separtion of variables
approach in which T(x,t)=F(x)*G(t). For the case where the
initial condition is T(x,0)=0 and the boundary conditions are
T(0,t)=1 and T(1,t)=0, one finds the solution shown in the
accompanying graph at the four indicated times. Note that t is
here non-dimensionalized via the ratio of the square of
the bar length divided by the thermal diffusivity.**

**TEMPERATURE
IN
A FINITE LENGTH BAR FOR DELTA FUNCTION INPUT: An interesting separation of variables
solution of the heat conduction equation occurs when one deals
with a finite length bar maintained at zero temperature at its
ends at x=0 and x=L and subjected to an initial delta function
d(x-L/2) input at its middle. The
solution for temperature in the bar will be worked in a
homework problem and found to have the value
T(x,t)=(2/L)*Sum[sin(n* p/2)*sin(n*p*x/L)*exp(-a(n* p/L) ^{2}*t),
n=0..Inf]. We show you here a MAPLE generated graph of this
solution for at=0.002 and at=0.02 , when L=1. As expected, this
solution looks a lot like the result for a bar of infinite
length when x is near L/2.**

**COOLING
OF
THE EARTH: An interesting
heat conduction problem concerns the cooling of a sphere of
radius r=a and an intial constant temperature T _{ o}.
This is a problem first looked at by Lord Kelvin in the 19th
century to determine the age of the earth. Casting the problem
into spherical coordinates one needs to solve T_{t}= a [T_{rr} +(2/r)T_{r}]
subject to T(0,t) finite, T(a,t)=0, and T(r,0)=T_{ o}.
Using the substitution T=R(r)/r]exp(- al^{
2}t), this leads to R"+l^{
2}R=0. From this follows the closed form solution
T(r,t)=Sum[(2T_{ o}a/npr)(-1)^{n+1}
sin(npr/a) exp(
a(np/a)^{2}t),
n=1..infinity]. We have plotted this result on the
accompanying graph for at=0.1, 0.5,
and 1. The approximate e-folding time for the original
temperature is seen from this result to be t*=(a/ p )^{2}/a.
Although Kelvin estimated from his solution (based on the
temperature rise in deep mines) that the earth was only some
24 million years old and this value is clearly in error as
pointed out by numerous sources at the time(Huxley etc), the
value of t*obtained for the earth ( assuming it to be made
essentially of iron where a=0.205cm^{
2}/sec) is actually t*=(6.378x10^{ 8}cm/ p )^{ 2}/0.205=2.01x10^{ 17}sec=6.38
billion years, which is in the right ballpark for the earth's
current estimated age of about 3.7 billion years and a bit
higher than this value because of the neglect of known
convection which speeds up the heat transfer process . Perhaps
Kelvin's shorter cooling time estimate was partially
influenced by the Victorian belief that , according to Bishop
Usher(1581-1656) , the earth was created precisely in 4004 BC.**

**TEMPERATURE RISE IN A BAR WITH SELF-HEATING: Another
problem of interest in the heat conduction area is that of the
temperature rise produced within a conductor when heat sources
are present within the conductor such as might be the
case for a nuclear reactor or a current heated wire.
Here the 1D heat conduction equation reads T _{t}=aT_{xx}+f(x,t) subject to the
BCs T(0,t)=T(L,t)=0, and an assumed IC of T(x,0)=0
. Here f(x,t) is the time and spatially dependent heat source.
Trying a standard sine solution of the form
T(x,t)=Sum[C(t)sin(n p x/L),
n=1..infinity] and also expanding f(x,t) as Sum[D(t)sin(n px/L), n=1..infinity)], leads to the
first order ODE dC(t)/dt+ a(np/L)^{2}C(t)=D(t). This
equation can be solved in closed form to yield the desired
solution. In the accompanying graph we show the result for the
case of a constant heat source f=1 at a=1
and L=1 when 30 terms are taken in the sum. The analytic form
of the solution reads T(x,t)=Sum[2*(L^{ 2}/a )*(1-(-1)^{n})*(1-exp(- a*(n p/L)^{2}t)*sin(npx/L)/(n p)^{ 3}, n=1..infinity]. Note
that as t approaches infinity, the steady-state solution is a
parabola and the above result without the exp term is just the
Fourier sine series for this parabola.**

**SOLUTION OF THE LAPLACE EQUATION IN THE
HALF-PLANE: We have
shown in class via a Fourier transform that the 2D Laplace
equation f _{ xx} +f_{yy}=0 in the half-plane
y>0, -inf.<x<+inf. , when subjected to the Dirichlet^{#
}boundary condition f(x,0)=f(x)
, is: f(x,y)=(y/ p)*Int[f(z)/(y^2+(x-z)^2),{-inf<z<inf}].
This integral has a particularly simple solution when f(x)
consists of delta functions. We show here the solution for
f(x)= d (x-1)+
d(x+1). You could view the resultant contours as the
isotherms under steady-state conditions for a semi-infinite
slab conductor subjected to two hot spots along the x axis at
x=1 and x=-1. Pierre-Simon Laplace(1749-1827) was perhaps
France's greatest mathematician ever. He was a professor at
the Ecole Militaire and later at the Ecole Normal, and also a
member of the French Academy of Sciences. His contributions
where numerous including work in potential theory, celestial
mechanics, probability theory, and Laplace
transforms. He was involved in the introduction of the
metric system and was appointed by Napoleon as Minister of
Interior. However, this last position lasted only a few
weeks when Napoleon realized that Laplace was not suited for
the job as "he was bringing the theory of infinitesimals into
the management of government ". Upon restoration of the
Bourbons, Laplace was made a marquis, but had lost most of the
support of his scientific colleagues because of his political
opportunism. You can read more about Laplace by going HERE
.**

**#-Johann Peter Gustav Lejeune
Dirichlet(1805-1859) was a mathematics professor at Breslau,
Berlin, and Goettingen. He is best known for the convergence
condition of a Fourier series for a function with
discontinuities.**

**SOLUTION
OF THE LAPLACE EQUATION IN A SQUARE: A very simple solution of the Laplace equation
occurs for the Dirichlet problem f _{
xx} +f_{yy}=0
subject to f(0,y)= f(x,0)= f(x,1)=0
and f(1,y)=1. A separation of
variables approach leads to the solution f
(x,y)=(2/p)*Sum[(1-(-1)^n)*sin(n py)*(sinh(npx)/(n*sinh(n p)), n=1..inf ]. We show you here the
contour plot for the iso-values of 0.75, 0.5, 0.25, 0.1 and
0.01 using a twenty term approximation to the infinite series.
Note, as expected from the theorem of the mean, the value of
f(0.5,0.5) at the square center is exactly 1/4=average value
of f on the four edges. You can
also solve the Laplace equation for more complicated boundary
conditions by the technique of superposition. For example,
with the Dirichlet conditions j(0,y)= j(1,y)=1 and j(x,0)=j(x,1)=0 you simply superimpose the
above solution and one where x is replaced by 1-x to get the
interesting harmonic function shown HERE.**

**SOLUTION
VIA FINITE FOURIER SINE TRANSFORMS: Certain problems require the solution of the
Laplace's equation in a semi-infite strip. Such problems are
conveniently solved by use of a finite Fourier sine transform
on the variable transverse to the strip axis. For the case of
f _{xx} +f_{yy}=0
in 0<x< p, 0<y<inf,
when the boundary conditions are taken
as f(x,0)=1 and f(0,y)=f(p ,y)=0, one finds that the finite sine
transform operation S[f(x)]=Int[f(x)*sin(nx), x=0..p] leads to an ODE solution
g(n)=(-1)^n*[exp(-n*y)-1]/n which is consistent with the
left boundary condition. Inverting this transform results in
the solution f(x,y)=(2/ p)*Sum[g(n)*sin(n*x) , n=1..inf],
which can be summed up to yield the function f(x,y)=x/p-(2/p)arctan[sin(x)/(exp(y)+cos(x))].
Iso-values
for this last function are shown in the plot.**

**FOURIER
SERIES SOLUTION OF THE LAPLACE EQ. IN A CIRCLE :We have shown that the solution to the 2D
Laplace equation within a circle r<a for a specified edge
condition T(a, q) can be
represented as a Fourier sine-cosine series multiplied by a
power term in (r/a)^n. For special cases of T(a, q ) this series will have only a few
non-vanishing Fourier coefficients and thus leads to very
simple analytic solutions. One such example occurs for T(1,q)=sin^2(q
)=0.5[(1-cos(2 q)] where the
solution is T(r, q)=(1/2)[1-r^2*cos(2 q)]=(1/2)*[1-x^2+y^2]. We have ploted
the contours for this function for the iso-values 0.1 through
0.9 at 0.1 intervals in the accompanying graph. Note that T at
the circle center and along two diagonal lines has the mean
value of [sin( q)]^2=1/2 as
expected from the theorem of the mean and the symmetry of the
problem.**

**SOLUTION
OF THE LAPLACE EQUATION INSIDE A CIRCLE VIA THE POISSON
INTEGRAL: -We have shown that a harmonic
function satisfying Dirichlet boundary conditions on a circle is
given by Poisson's Integral. The integral reads
Phi(r,theta)=Int[(a^2-r^2)*f(u)/(a^2-2*a*r*cos(theta-u)+r^2),{u,0,2*Pi}],
where
a is the circle radius, u the integration variable, and
f(u) gives the value of Phi on the circle at r=a. The explicit
solution for Phi, for the special case of a delta function
Dirichlet condition , is shown in the accompanying graph.**

**FLOW ABOUT A CORNER: Consider
next inviscid flow about a corner of angle q. Solving the
Laplace equation for the steamfunction Psi and demanding that
Psi vanishes on the walls intersecting at the corner, one finds
the exact solution Psi=r^b*sin(b*Theta) with b=Pi/q. If we now
take the special case of a corner flow with q=Pi/4 radians, this
leads to Psi=r^4*sin(4*Theta) which in Cartesian coordinates
reads Psi=4*x*y*(x^4-y^4).**

**A contourplot ( looking very much like a spider web) is shown.**

**STREAMFUNCTION
AND
VELOCITY POTENTIAL FOR A DOUBLET: -We have shown
that 2D irrotational flows can be represented in terms of a
complex velocity potential F(z)=f + i
y, where fis
the velocity potential and y is the
streamfunction of the flow. We show here the pattern for the
doublet F(z)=1/z which is equivalent to f=x/(x^2+y^2)
and y =-y/(x^2+y^2). Note that f=const. and y =
const. always forms an orthogonal familiy of curves and are each
a solution of a 2D Laplace equation. This follows directly from
the Cauchy-Riemann conditions for any analytic function.**

**HARMONIC FUNCTION FOR INVISCID FLOW ABOUT A
CYLINDER: -The streamfunction for irrotational
flow about a cylinder satisfies the Laplace equation plus the
boundary conditions that the normal velocity at the cylinder
surface vanishes and that at large distance from the
cylinder the velocity is constant and parallel to the
x-axis. The solution of this boundary value problem for the
streamfunction y (x,y) is shown in
the attached and can be constructed by simply adding a doublet
(1/z) to a rectilinear flow (z). The heavy blue curve represents
the streamfunction value of zero located on the cylinder and
along the x axis. The stagnation points occur were the
x-axis intersect the cylinder.**

**SPIRAL FLOW: Another solution to
the 2D Laplace equation for inviscid flow is the harmonic
function Psi(x,y)=0.5*ln(x^2+y^2)-m*arctan(y/x). This solution
represents the steamlines Psi(x,y) for a spiral flow whose
complex velocity potential can be expressed as
F(z)=(-m+i)*ln(z). Here m is the arctan of the spiral pitch
angle. When m=0 one has a simple vortex flow while for m equal
to infinity one has a source flow. Some interesting spiral flows
encountered in nature are shown HERE and HERE.**

**2D FLOW PRODUCED BY SUPERIMPOSING TWO F(z)s:Another solution of the Laplace equation is
produced by superimposing the rectilinear flow solution F _{1}(z)=U_{o}z
and a source flow solution F_{2}(z)=[Q/(2p)]ln(z). in this case the combined
streamfunction is y(x,y)=U_{o}y+[G/(2p)] arctan(y/x). A contour graph for
these streamlines when G/(2pU_{o})=1 is shown here and is
seen to corresppond to flow about a blunt body. One may also
construct an infinite number of other combinations. For
example , a source and a sink at z=+i and =i, respectively,
plus a superimposed counterrotating vortex at z=1 and
co-rotating vortex at z=-1 leads to the streamline pattern
shown HERE when
Q=G.**

**INVISCID FLOW
ABOUT A ROTATING CYLINDER:An
interesting extention of flow about a stationary cylinder
occurs when one superimposes a vortex flow to get the complex
velocity potential F(z)=U _{o}(z+a^{2}/z)-i( G/2p )ln(z). In this case the front
and rear stagnation points are moved away from the x axis and
the cylinder will experince a lift at right angles to the flow
direction. We show you here the velocity streamfunction
pattern y(x,y)=Re[F(z)]=Constant for a special case. Note the
high velocity region above the cyliner(ie. close streamline
spacing) and the lower velocity region below the cylinder.
From Bernoulli we know that this will result in a lift( Magnus
Effect). It is interesting to note that this solution can be
conformally mapped via a Joukowsky mapping to that of flow
about an airfoil and hence predict the latter's lift(see
below).**

**STREAMLINE
PATTERN
PRODUCED BY TWO COUNTER-ROTATING VORTEXES: Another interesting inviscid flow pattern can
be constructed by superimposing the effect of two vortexes. We
know from our class discussion that a single vortex of
circulation G placed at z _{o}=a+ib
in the z plane, has its streamfunction given by y(x,y)=-(G/4p)ln[(x-a)^2+(y-b)^2)].
Consider
now two vortexes one placed along the x axis at a=1 and a
second one of opposite circulation placed along the x axis at
a=-1. In this case one readily finds that y
(x,y)=( G/2p)ln{[(x+1)^2+y^2]/[x-1)^2+y^2]}.
As the plot shows, this streamline pattern for y (x,y)=Const consists of a collection
of circles centered along the x axis. Note the strong downward
velocity component V=- y (x,y)_{x
}along the y axis. Such patterns are very reminiscent of
the wind patterns observed between a high and a low pressure
region on a weather map. In the northern hemisphere the
circulation direction about a high is clockwise and accounts
for the fact that when we do get a few days of cold weather in
the wintertime here in Gainesville it usually means that a
high pressure region is located directly north of us ,
such as over Columbus,Ohio.**

**CONFORMAL
MAPPINGS OF 2D REGIONS: An interesting property of the 2D Laplace
equation is that it is invariant under a conformal mapping
w=u+i*v=f(z). Thus one is able to solve it in quite
complicated regions by transforming this region into a
simpler one where an analytic solution is known and then
translating back. As an example of such, consider trying to
solve the Lapace equation inside a cardioid whose boundary is
given by R=2*[1+cos( Q)] . Here
R=sqrt(u^2+v^2) and tanQ =v/u. This
would be clearly difficult to solve in its present
configuration, but becomes quite easy to treat after applying
the mapping w=R*exp(i* Q)=z^2 ,
where z=x+i*y. In this case the image of the cardiod becomes
the shifted circle of unit radius r=2*cos(
q ) and the problem can be solved by the Poisson
integral solution discussed earlier. We show you here the
cardiod and its circle image in going from the w to the z
plane.**

**JOUKOWSKI AIRFOILS:- The
figure shows the airfoil like shapes obtained by applying the
Joukowskimapping f(z)=z+1/z to an off-center circle
defined by (x+a)^2+(y-b)^2=[1+sqrt(a^2+b^2)]^2. This
transformation allowed the earliest calculation of airfoil lift
by reducing the problem to one inviscid flow about a rotating
cylinder. As is seen, the form of the closed curves produced by
the transformation depends very much on where in the second
quadrant the circle center (x=-a, y=b) is located. A very nice
applet allowing you to see the image produced by using different
circle centers is found by clicking HERE.
Nikolai Joukowski (1847-1921) was a professor of mechanics
and applied mathematics at the Moscow Technical School and
received his PhD from Moscow University. He wrote some 200
papers in his lifetime and is considered the father of the
Russian School of hydrodynamics and aeromechanics.(I thank
Raymond Deleu of the Netherlands for supplying me with the
information on Joukowski. The reason I did not have this
information before is that I was not using the
correct Russian spelling "Zhukovskii" when searching.)**

**ISOTHERMS
IN
AN INFINITE STRIP: A very nice application of
conformal mapping is the use of w=exp(z) to determine the
steady-state temperature in an infinite strip
-infinity<x<+infinity, 0<y<p when T(x,y) vanishes everywhere
on the boundary except for that part between x=-1 and x=+1
where T(x,0)=1. This is clearly a difficult problem to solve
in the x-y plane but becomes quite manageable when converting
to the w=u+iv plane. Using the conformal mapping indicated,
one finds u=cos(y)*exp(x) and v=sin(y)*exp(x). In the w plane
the strip maps into the upper half plane u>0 and one can
now use the earlier derived integral solution of the Laplace
equation. Doing so, this yields T(u,v)=(1/p)*[arctan((e-u)/v)-arctan((1/e-u)/v)]
.
Substituting in for u(x,y) and v(x,y) one gets the result and
isotherms shown in the attached graph.**

**SOLUTION
OF
THE LAPLACE EQUATION IN A SEMI-INFINITE STRIP: A problem which we have considered earlier
using finite Fourier transforms may also be solved by a
conformal mapping. The problem is finding the harmonic
function in the semi-infinite strip 0<x<infinity,
0<y<p which satisfies the Dirichlet boundary
conditions f(0,y)=1, f(x,0)=0
and f(x, p)=0. The solve this problem we
use the conformal mapping w=cosh(z) so that
u=cos(y)cosh(x) and v=sin(y)sinh(x). This mapping converts
the semi-infinite strip into the upper half of the z plane.
there it is a simple matter to use the known Poisson
integral solution for the half-plane derived earlier.
integrating this integral yields, after a little
manipulation , the final result that the harmonic function
is f=(2/p )arctan(sin(y)/sinh(x)]. A
contour plot of this function is shown in the accompanying
graph for the values 0.9, 0.8, 0.7 etc.**

**APPLICATION
OF THE SCHWARZ-CHRISTOFFEL TRANSFORMATION: The
Schwarz-Cristoffel transformation is a particular type of
complex variable mapping which opens up an n sided polygon
within the w=u+i*v plane into a half-plane within the z=x+i*y
plane. This allows one to obtain exact analytic solutions of the
Laplace equation within the polygon for an imposed set of
Dirichlet boundary conditions by using the known solution of the
Laplace equation in the half-plane given by the Poisson
integral. We demonstrate this procedure for a solution within
the sector 0<R<infinity, 0<Y<p/4
given the bcs that T=0 along the boundary except for that part
from 0 to +1 along the lower boundary where T=1. For this case
the SC transform follows from w=A*Int[z-0)^(-3/4)]+B, which
yields w=z^0.25 , after requiring that the image point of w=0 is
z=0 and of w=1 is z=1. The Laplace equation in the z-plane now
has the known solution T(x,y)=(y/p)*Int[dx
/(y^2+(x-x)^2),{x,0,1}]which integrates to
T(x,y)==(1/p)*{arctan[(1-x)/y]+arctan(x/y)}. Using the double
angle formula for arctan, this leads to
T(x,y)=(1/p)*arctan[y/(y^2+x*(x-1))], from which the desired
solution T(u,v) within the sector follows via the equality
(u+i*v)^4=(x+i*y). A drawback of the SC transform is that
in most instances one cannot obtain an analytic solution to the
SC integral representing w(z) . For those cases were exact
integrations are possible, like the one presented in the present
example, we refer you to the Dictionary of Conformal
Representations by Kober. Karl Schwarz(1843-1921) was a
professor at Halle, Zurich, Goettingen and Berlin specializing
in the calculus of variations and conformal mappings. His name
is also attached to the Schwarz inequality. Elwin
Christoffel(1829-1900) was a mathematics professor at Berlin,
Zurich, and Strasbourg working on function theory, tensor
methods and differential equations.**

**CONFORMAL
MAPPING TECHNIQUES USED TO DETERMINE THE STREAMFUNCTIONS FOR
INVISCID FLOW PAST A PLATE BARRIER: We have shown in class that the Laplace
equation is invariant under the transformation w=f(z) in going
from the z=x+i*y to the w=u+i*v plane . Furthermore the angle
is conserved (hence the transformation is conformal) as long
as df(z)/dz does not vanish. The
Schwarz-Christoffel transformation is a particular complex
function mapping which opens up an n sided polygon into a half
plane. We can use it to solve, for example, the problem
of inviscid flow past a plate barrier. In this case the
rectilinear flow in the half plane z>0, for which the
streamfunction is y=y, can be
transformed into a barrier plate configuration using the SC
mapping w=sqrt(z^2-1). It leads to the streamline pattern
shown in the figure. Analytically this flow field can be
expressed as y^4+ y^2*(1+u^2-v^2)-(u*v)^2=0 or in its
equivalent form y=sqrt[0.5*(-1-u^2+v^2)+sqrt(0.25*(1+u^2-v^2)^2+(u*v)^2)]
,
where u and v are the real and imaginary parts of w. Note that
for a real fluid, where viscosity is present , separation will
occur at the plate edges so that a wake producing drag is
established on the downstream side of the flow. For the pure
inviscid flow considered here the drag will be zero.**

**SPHERICAL
HARMONICS:In solving the 3D Laplace equation in
spherical coordinates(r, q, f) one
finds that the angle portion of the solution is given by the
spherical harmonics Y(q,f
)=Const.*P[n,m](cos(q))*exp(i*m* f). Here P[n,m] is the associated
Legendre polynomial related to the normal Legendre
polynomials by P[n,m]=(1-x^2)^(m/2)*d^m[P[n](x)]/dx^m with
x=cos(q ) and P[n] given by the
Rodrigues formula P[n]=1/(2^n*n!)*d[(x^2-1)^n]/dx^n. We
treat q as the polar angle and f as the azimuthal angle. This is the
standard notation used by physicists and engineers , but is the
reverse of what is used by mathematicians and can thus lead to
confusion. We look here at one particular example, namely the
spherical harmonic Y(4,3). To construct this function our
starting point is P[4]=(1/8)(35*x^4 -30*x^2 +3) which, after a
simple manipulation, leads to P[4,3]=150*x*(1-x^2)^1.5= 150*sin( q)^3*cos( q
). Thus, taking the real part, one finds that Y[4,3]=Const.*sin( q)^3*cos( q)*cos(3*f). The graph shows a 3D plot of this
function expressed as r=abs[Y(4,3)] obtained by using sphereplot
in MAPLE. Since any spherical surface function F(q,f) can be expressed as an infinite sum
of spherical harmonics, the Y's are encountered in a variety of
applications including geodesy , quantum mechanics, and quantum
chemistry.These harmonics are also being used by computer
artists to construct some interesting looking 3d images. By
clicking HERE
you can see some thousand different spherical harmonic
shapes, or look at something I drew using Y[8,4] and
windows paintbrush by clicking HERE .**

**STEADY
TEMPERATURE IN A SPHERE EXPRESSED VIA SPHERICAL HARMONICS:
We have shown in class that the
solution of the Laplace equation inside a sphere of radius r=a
is T(r,q,j)=Sum[A _{nm}(r/a)^{n}*Y_{n}^{m}(q,j), n=0..infinity, m=0..infinity].
This solution becomes particularly simple when there is no
azimuthal angle dependence of the surface condition at r=a.
Under that limit T(a,q,j)=f(q) which means m=0 and your spherical
harmonic term reduces to the standard Legendre polynomial P_{n}(cosq). As an example take the case of the
steady-state temperature distribution within a sphere of
radius r=a=1 when f(q)=1-(cosq)^{2 }=0.5(1-cos(2q). Here there are only two
non-vanishing terms in the infinite double series and these
non-vanishing coefficients, determined with aid of the
orthogonality relation Int[P_{n}(x)P_{m}(x),
x=-1..1]=2/(2n+1), are A_{0,0}= 2/3 and**

**GREEN'S
FUNCTION FOR THE SOLUTION OF THE LAPLACE EQUATION WITHIN
THE HALF-SPACE:-We have derived in class the fundamental
formula for the solution of the Laplace equation in a 3D space.
By use of the Kelvin charge method one is able to construct
Green's and Neumann's functions which allow simple solutions for
certain types of 3D geometries. One of these problems deals with
the solution of Laplace's equation in the half space z>0,
given the Dirichlet boundary condition f(x,y,0)=f(x,y).
Here the appropriate Green's function can be constructed by two
point charges placed symmetrically above and below the x-y plane
and having opposite unit charge. That is G=(1/4p)*[1/r-1/r') where r and r' are the two
distances from the point charges to the movable point
Q(x,y,z) . If we look at the cross-section of this function
found by cutting it with the y=0 plane, one finds the graph
shown. It clearly indicates that G vanishes everywhere on the
x-y plane as required. George Green(1793-1841) was a
mathematical genius whose formal education had stopped with
grade shool. His occupation as a mill owner left him sufficient
time to write ten or so remarkable scientific papers. His paper,
"An Essay on the Application of Mathematical Analyisis to
Theories of Electricity and Magnetism" , brought him to the
attention of the scientific establishment and he was invited,
and then began, to attend Cambridge University as an
undergraduate at the ripe old age of forty. Unfortunately his
health deteriorated and he died a few years later . His name is
associated with Green's theorem , the various Green's formulas,
and he is also credited with invention of the Green's function.
He was never married but had seven children with the daughter of
his mill foreman.**

**NEUMANN
FUNCTION FOR THE INFINITE HALF SPACE PROBLEM: -The fundamental formula for the solution of
the Laplace equation in 3D can also be used as a starting
point for solving a Neumann boundary value problem. Consider
again the half space z>0 and this time impose the boundary
condition that d f (x,y,0)/dn
=g(x,y) is given. This time the Kelvin charge method requires
two equal sign charges, one placed above and the other
symmetrically below the x-y plane. The resultant Neumann
function reads N=(1/4 p)*[1/r+1/r']
and the Laplace equation solution becomes
f(x,y,z)=Int(N*g(h,x)*d hdx). In
the accompanying graph we show the contour map of the
intersection of the y=0 plane and the N=Const surfaces .
Note the zero values of dN/dn at all points on the x-y plane.
These projections for N=Const. are reminiscent of the Ovals
of
Cassini defined slightly differently by
r*r'=Const. Karl Gottfried Neumann(1832-1925) was a professor
at Leipzig, the same university at which my father got his
physics PhD in 1933 and where Heisenberg(uncertainty
principle) and Teller(future H bomb builder) were his
contemporaries. HERE is a photo of K.G. Neumann in his later years.
Looks like he needs a haircut.**

**SOLUTION OF
THE LAPLACE EQUATION IN THE HALF-PLANE USING A GREEN'S
FUNCTION:- Earlier
we had shown that one can use a Fourier transform approach to
solve the Laplace equation in the half-plane y>0 knowing the
value of the harmonic function everywhere on the x axis. We show
you here a second way to arrive at an identical result
using the Green's function route. The fundamental theorem in 2D
allows one to write f (P)=--Int[ f(C) dG/dn dl], where C is the border of
a closed region and P is any point within this region. Here dl
is the increment of length along the bounding curve C and dG/dn
the partial derivative of a Green's function with respect to the
outward normal to the curve C. For the special case of the
half-plane , C becomes the x axis , P is located at ( x _{0},y_{
0}) in y>0, and G is constructed from two line charges
of opposite sign placed at P(x_{0} ,y_{ 0)}and
P'(x_{0},y_{ 0}). The explicit form of Green's
function becomes G=-(1/(2*Pi))*[ ln(r)-ln(r')], where r and r'
are the distances from P and P' to a point Q on the x
axis. The explicit form of the Green's function plotted in the
graph for x_{0}=0, y_{ 0} =1 is
G=-[1/4*Pi)]*ln{[(x-x_{ 0})^2+(y-y_{ 0})^2]/[(x-x_{
0})^2+(y+y_{0} )^2]}. The solution of the Laplace
equation thus becomes the Poisson integral f
(x_{ 0},y_{0})=[(1/Pi)]*Int[ f(x) dx/((x-x_{0}
)^2+y_{ 0}^2), x= -inf..inf], where f
(x) is the x dependent value of f along
the
x
axis.**

**SPERICAL BESSEL FUNCTIONS OF THE FIRST KIND
[j(n,x)]: - These functions arise in a solution of the
Helmholtz equation in spherical coordinates and are governed
by the ODE x^2y"+2xy'+[x^2-n(n+1)]y=0. One has that
j(n,x)=sqrt[Pi/(2x)]*J(n+1/2,x) so that j(0,x)=sin(x)/x.
Hermann von Helmholtz (1821-1894) was a professor of anatomy
and physiology at Bonn and Heidelberg and later professor of
physics at Berlin. His name is also associated with the
Kelvin-Helmholtz instability in fluid mechanics, the
conservation of energy law, an acoustic resonator, the
ophthalmoscope, and the theory of hearing. You can find a
picture of HvH by going HERE.**

**SPERICAL
BESSEL FUNCTIONS OF THE SECOND KIND[n(n,x)]: -The second set of functions satisfying the
spherical Bessel equation given above. These solutions all
approch minus infinity as x-goes to zero. This is to be
expected since x=0 is a regular singular point of the
governing equation. The identity
n(n,x)=(-1)^(n+1)*sqrt[Pi/(2x)]*J(-n-1/2,x) holds. From this
follows that n(0,x)=-cos(x)/x.**

**TRANSMISSION
AND
REFLECTION OF A SOUND WAVE AT AN INTERFACE: -A one dimensional acoustic wave approaching
an interface will be partially reflected and partially
transmitted . The governing equation for the case of a
strictly periodic wave is the Helmholtz equation y _{xx}+k^{2} y=0 with solutions y_{R}=A_{ R}exp(ik_{1}x)
and y_{I} =A_{ I}exp(-ik_{1}x)
for the reflected and incoming wave at the left of the
interface. Also the solution is y_{
T} =A_{T}exp(-ik_{2}x) for the
transmitted wave to the right of the interface. For all three
waves the time dependence goes as exp(i wt).
Now at the interface both the velocity of the fluid elements,
y_{x }, and the
pressure , which by Bernoulli reads -i
rwy, are continuous for the small amplitude waves
considered. Eliminating A_{T} from the two
resultant equations, yields the reflection coeficient R=(A_{T}
/A_{I})^2 , which, after a short manipulation, can be
expressed as R=[(1-g )/(1+ g)]^2 and represents the fraction of
the incoming wave energy which is reflected. Here g =( r_{
2} c_{2})/( r_{ 1}c_{
1} ) is the acoustic impedance ratio with r being the density and c the speed of
sound in the same material. The transmission coefficient is
that fraction of the wave energy penetrating the interface and
is simply equal to T=1-R=4 g /(1+ g )^2. The plots for R and T are shown
in the graph and indicate large reflections whenever g differs a lot from unity. At a
water-air interface one has g=1.5*10^5/42.8=3504.67
meaning
some 99.9% of the signal is reflected. It is this poor
transmission of acoustic signals across a water-air interface
which makes it difficult to pick up sonar signals from a
submerged submarine from an airplane without the presence
of ocean surface transponders.**

**SOUND TRANSMISSION THROUGH A BARRIER:An extention of the Helmholtz equation
solutions for sound waves at an interface is concerned
with the transmission of a sound wave through a finite width
barrier. This time one has a bit more complicated problem
involving velocity and pressure continuity at two interfaces.
Click HERE
to see a full development for
the transmission coefficient T=1-R across a layer of thickness
d. Note that the graph shown applies to the special case of
transmission through a layer of glass bounded on both sides by
air. The results are quite interesting showing very poor sound
transmission as long as the acoustic impedance ratio is far
removed from unity. Also one observes from our equation
the known contructive interference phenomenon when the sound
wavelength in the glass equals multiples of the barrier
thickness. These points however are difficult to reach unless
the product of the barrier thickness and sound frequency is
comparable to the sound speed in the barrier.**

**FRAUENHOFER
DIFFRACTION
BY A SLIT: We have shown
in class that the fundamental formula for the solution of the
Helmholtz equation when applied to a light path from a souce S
passing through a slit and onto a photographic plate P is
given by the wavefunction f(P)=Const.*Int[exp-i*k*(r+R)
dS , where k is the light wavenumber, dS an increment of the
slit aperture, r the distance from a point in the slit to the
plate P, and R the distance from the light source S to the
same point in the aperture. This result is valid as long as
the light wavelength is short compared to the aperture(slit)
width of 2a. Expanding r+R in a Taylor series about the
aperture center yields the various diffraction orders from
first order(Frauenhofer) through second order(Fresnel) etc.
For a linear slit, the lowest order diffraction yields f(P)=Const.*Int[exp(-ikd qh),h=-a..a],
where dq is the difference in angle
between the incoming and the diffracted ray. Solving this
integral yields the Frauenhofer diffraction result f(P)=Const *(sin(x)/x), where x=k*d q*a.
We show this pattern in the graph. Note that one cannot
resolve two closely spaced sources if the subtended angle
between them is closer than x=p
or, equivalently, d q =p/(k*a)=l
/2a. One is said to have reached the diffraction limit of an
optical instrument when this angle is reached. It explains one
of the reasons why good astronomical telescopes use large
aperture lenses or mirrors. Also why an electron microscope,
where the equivalent wavelength l can
be as low as several Angstroms as compared to about one-half
micron for visible light, has a much better resolution than an
optical one.**

**GRAVITATIONAL
POTENTIAL
INSIDE AND OUTSIDE A UNIFORM DENSITY SPHERE:- It is known that the gravitational potential f is determined by solving the Poisson
equation dell^2 f=-4*Pi*G*r, where G is the universal
gravitational constant and r the
spatially dependent density r(x,y,z).
When dealing with an infinite space, the Green's second
theorem allows one to obtain the closed form solution f(x _{o},y_{o},z_{ o})=G*Int[ r(x,y,z)dxdydz/r], where r=sqrt[(x-x_{o})^2+(y-y_{
o} )^2+(z-z_{o})^2] and the integration extends
over all values of x, y and z. If we now consider a uniform
density sphere of radius r=a, then the above integral can be
solved exactly to yield f =M*G/z_{o}
for any point on the z axis at distance z_{ o}>a
away from the sphere center( Newton's result). Here M is the
sphere mass. Inside the sphere one has to break the integral
up into two parts, one for 0<r<z_{ o} and the
second for z_{ o}<r<a. A short manipulation then
leads to**

**LEGENDRE
MULTIPOLE EXPANSION: For all except
the simplest geometries, the integral representing the
gravitational potential solution V(P) to
(dell^2)V=-4*Pi*G*r cannot be
evaluated in closed form. Under this more general condition,
one makes use of a multipole expansion method first introduced
by Legendre. The idea is to expand the term T=1/sqrt[R _{o}
^2+r^2-2*r*R_{ o}*cos(q)]
as T=(1/R_{ o})*[1-(1/2)* e+(3/8)* e^2-...], where e=(r/R_{o}
)*[(r/R_{ o})-2*cos( q )].
This will be a rapidly converging series whenever e<<1 and thus when the distance
R_{o }from the observation point to the body's mass
center is large compared to all distances r from the
mass center to any point in the mass(see figure).
On expanding this result in powers of r/R_{o} ,
one finds the solution V(P)=(G/R_{ o})*Int[r*Sum((r/R_{ o} )^n*P(n, cos(q), n=0..inf)dV], where P(n,x) are the
Legendre polynomials. As you probably recall, they have the
values P(0,x)=1, P(1,x)=x, P(2,x)=(1/2)*(3*x^2-1) etc. The
first term in this series approximation for V(P) represents
the monopole GM/R_{o} , which will be the only
non-vanishing term for a spherically symmetric mass
distribution such as a self-gravitating , but non-rotating ,
star. The sun is believed to also have an additional small
quadrupole (n=2) term due to a mass asymmetry produced by its
rotation.**

**THE LAGRANGE POINTS FOR THE EARTH-MOON
SYSTEM: We have developed
the general solution for the gravitational potential of a
distributed mass. Let us now use this result to look at the
famous problem of Lagrange(1772) dealing with the potential
field produced at a point P in the neighborhood of two bodies
rotating about their common center of mass. We'll consider the
earth-moon system where M=5.97x10^24 kg and m=7.349x10^22 kg
so that M/m=81.235. Letting the distance from the earth center
to the cg point be x _{1} and that from the cg to the
moon center be x_{2}, we choose a length scale where x_{
1} +x_{2}=1 for the earth-moon distance. Now
the total gravitational potential felt at a point P at
distance r_{o} from the cg is GM/r_{M}+Gm/r_{
m} , where r_{ M }is the M-P distance and r_{m }
is the m-P distance. We also know from a balance between
Newton's law of attraction and either the earth's
centrigugal force M w^{2}x_{
1} or the moon's centrifugal force mw^{2}x_{ 2 },
that w^{2}=GM/(x_{2}(x_{
1} +x_{ 2})^{2})**

**j(x,y)=GM{1/sqrt[x^2+y^2+x _{
1}^2+2*x_{1*}x]+(m/M)/sqrt[x^2+y^2+x_{2}^2-2*x_{2}
*y]+(x^2+y^2)/(2*x_{2}*(x_{1}+x_{2}
)^2)}.**

**A contour plot of this function for the
earth-moon case is given here. As seen the points L1, L2, and
L3, lying along the earth-moon axis, exist at saddle points of
the contour plot and hence are unstable, while the L4 and L5
points lie at a potential minimum and thus are stable. These
stable points occur near j(x,y)/GM=1.5135
and they are located at the vertex of an equilateral triangle
having the earth and the moon located at the other two
vertexes. Until its
merger with the National Space Society(NSS) several years ago,
there existed a L5 Society dedicated to the future of man in
space. The L5 designation came from this fifth Lagrange point
which exists for any rotating double mass system when
M>>m. L5 would be a good place to assemble equipment for
a planetary trip since the gradient of the potential, and
hence the force, is zero there.**

**CONSTRUCTION OF BIHARMONIC FUNCTIONS:The biharmonic equation arises in a variety of
areas including elasticity and creeping flow. In the latter
case the governing equation for a 2D incompressible flow
reads y _{xxxx}+2***

**DISPERSION
LAW
FOR LINEAR WATER WAVES: For
the benefit of our coastal engineerinmg students, we
consider near the end of our PDE course the problem the
properties of water waves. One considers small amplitude
sinusoidal water waves of frequency w
and wavenumber k on the surface of a water layer of depth h.
Starting with the wave equation for an irrotational fluid ,
one can assume a solution for the velocity potential of the
form f=Real[F(y)*expi( wt-kx)], where x is the direction of
wave propagation and y the coordinate normal to the water-air
interface. Also the wave surface can be described by h=h _{ o}*expi( wt-kx). Substituting f into the wave equation f_{tt} =c^{2}( f_{xx}+ f_{yy})
, assuming c to be infinite (actually c is about 1500
meters per sec in water), and applying the bottom boundary
condition that the normal derivative of f
must vanish at y=-h, one finds the solution
F(y)=A*cosh(k*(y+h)). At the water-air interface , one
requires that both the velocity component in the y direction
and the pressure be continous there. That is d f/dy=d h/dt
and d f /dt=g*
h at h=0, with the
second equality following from the Bernoulli law for small
amplitude waves. We thus find that h_{
o}i w =Ak*sinh(kh) and i wA*cosh(kh)=-g h_{
o} , from which one obtains at once the well known
dispersion relation**

**THE
SCHROEDINGER EQUATION AND THE QUANTUM WELL PROBLEM: In 1926, E.Schroedinger converted the
classical conservation of energy equation E=p ^{2}/(2m)+V
to quantum mechanical language by introducing a wavefunction
having the form Y=yexp(-i
E t /hbar). Using the transformations ( based to some extent
upon experimental observations) that the momentum operator
p=mv becomes -i hbar dell(Y), the
energy transformation E->i hbar Y_{t},
and the potential transfomation V->VY,
he arrived at the
time-dependent Schroedinger equation i hbar Y_{t}=-(hbar)^{2}/(2m)*(dell^{2}Y)+VY. Here hbar is Planck's constant divided by 2p. Upon using the above stated form for
the wavefunction, this equation simplifies to its familiar
time-independent form, which in 1D, reads, y"+2m/(hbar)^{2}[E-V]y=0. Solutions to this equation
subjected to boundary conditions that the wavefunction
should vanish at infinity become possible only for certain
quantized energy levels E and the probability of finding an
electron of mass m at a given position x in a
potential field V(x) is proportional to the square of y.
To demonstrate the use of this equation, consider the
simple example of an electron sitting in a 1D potential well
- a<x<a where V=0. Outside the well the potential is
+V_{o }. In this case the even solution of
the Schroedinger equation in the well has the standing wave
form y=A*cos[sqrt(2mE)x/hbar] . Outside the well
for x>a the solution is y=B*exp[-sqrt(2m(V_{o}-E))x/hbar].
If
one now demands continuity of y
and its derivative
at x=a, one arrives at the transcendental relation sqrt[(V_{o}-E)/E]=tan[sqrt(2mEa^{2}/hbar^{2})]. From the plot of this result shown, one
sees that the energy levels associated with the electron in
the well can only take on certain quantized values. For
example, when the well becomes infinitely deep , so that V_{o}
approaches infinity, the lowest energy level will be at E= p^{2 }hbar^{2}/(8ma^{2}).
By clicking HERE you
can see the even wavefunction for the quatum well problem
corresponding to the third energy level when
sqrt(2mVo)*a/hbar=10, a=1 so that sqrt(2mE)*a/hbar=7.08.
Erwin Schroedinger (1887-1961) was an Austrian physicist who
held positions at Zurich, Berlin, Oxford, Graz, and the
Institute for Advanced Studies in Dublin. He won the Nobel
Prize for his wave equation work in 1933, served on the
Italian and Hungarian fronts during WWI, had a rather wild
personal life involving having children with at least three
different women while married to his original wife Anny, and
went somewhat off of the deep end later in life claiming to
have discovered a unified field theory.(So far no-one,
including major attempts by Einstein and Heisenberg, has
been able to come up with such a theory combining nuclear,
atomic, electromagnetic, and gravitational forces. It
remains the elusive Holy Grail for physicists)**

**EGM 6323-INTEGRAL
EQUATIONS AND VARIATIONAL METHODS**

*Integral equations of Volterra and Fredholm. Inversion of
self-adjoint operators via Green's functions. Hilbert-Schmidt
theory and the bilinear formula.The calculus of
variations. Geodesics, the Euler-Lagrange* *equation, and
the brachistochrone. Variational treatment of Sturm-Louiville
problems. Fermat's principle.*

**VOLTERRA EQUATION SOLUTION VIA NEUMANN
SERIES:-We have shown by
direct integration that a differential equation subjected to
initial conditions can be re-written as an integral equation
of the Volterra variety and that this type of integral
equation is solved by a series expansion of the Neumann type.
We here demonstrate the use of such a series in solving the
integral equation y(x)=x-x*Int[y(t)/t,{t,0,x}] starting with
the initial term y0=f(x)=x. It then follows that
y1=y0-x*Int[1,{t,0,x}]=x-x^2 and
y2=y1+x*Int[t,[t,0,x}]=x-x^2+x^3/2. It is clear that this
series converges toward y(x)=x*exp(-x) which also satisfies
the original differential expression y'+y=exp(-x) with the ICs
y(0)=0 and y'(0)=1 from which the IE arose. Karl Gottfried
Neumann(1832-1925)was a professor at Leipzig specializing in
mathematical physics, potential theory and electrodynamics.
The Bessel function of the second kind Y(n,x) is sometimes
referred to as the Neumann function N(n,x).**

**SOLUTION OF
VOLTERRA EQUATIONS WITH DIFFERENCE KERNELS :Volterra
equations of the type y(x)=f(x)+ l*Int[K(x-t)*y(t),t=0..x]
can
be solved exactly by use of Laplace transforms and application
of the Laplace transform of a convolution integral. A straight
forward procedure leads to y(x)=L ^{-1} [fbar/(1-l*Kbar)],
where the bar designates the L transform of the indicated
function. We consider here the special case of f(x)=x^2 , l=1, and K(x-t)=x-t. It leads to**

**SOLUTION OF A VOLTERRA EQUATION BY THREE
DIFFERENT METHODS: We have
shown in class that a Volterra equation can be solved by
(a)Neumann series, (b)Solution of the corresponding
differential equation, and(c) by Laplace transform methods .
The last technique works only for difference kernels. Consider
now the integral equation y(x)=x+Int[sin(x-t)*y(t), t=0..x].
This equation can be solved by each of the mentioned methods.
The quickest way to its solution is to use Laplace transforms
which yields ybar=1/s^2+1/s^4 and on inverting produces the
exact solution y(x)=x+(1/6)*x^3. The differential equation
route, done by differentiation and use of the Leibnitz rule,
produces the intial value problem y"=x, y(0)=0 and y'(0)=1 .
This yields the same exact solution. Finally the Neumann
series approach yields (with help of MAPLE) the series-**

**y(x)=x+[x-sin(x)]+[x+x*cos(x)/2-(3/2)*sin(x)]+[x-(15/8)*sin(x)+(7/8)*x*cos(x)+x^2*sin(x)/8]+**

**Click on the above title to see a plot of
the exact solution compared to the four term Neumann series
approximation. Vito Volterra(1860-1940) was a professor of
mechanics and later mathematics at the universities of
Pisa, Turin, and last at Rome. He lost his job in the early
1930s during the Mussolini regime. Click HERE
to see a photo of Volterra as called to my attention by Thomas
Chestnut.**

**TRIANGULAR KERNEL:-This is the
two variable symmetric kernel T(x,t)=x(1-t) for x<t and
T(x,t)=t(1-x) for x>t and folows from d^2y/dx^2=0 subject to
boundary conditions y(0)=y(1)=0, to the Fredholm integral
equation y(x)=Int[T(x,t)f(y,t),{t,0,1}]. For a derivation of
T(x,t) click HERE.Erik
Fredholm(1866-1927)
was a professor of mathematical physics at Stockholm and is
mainly known for his contributions to integral equations and
spectral theory. He wrote relatively few , but well received,
papers during his lifetime. Click HERE
to see a picture of Fredholm.**

**SOLUTION
OF THE FREDHOLM EQUATION VIA PICARD ITERATION: -A very powerful technique for solving the
Fredholm integral equation, especially when dealing with two
part Greens functions as the kernel, is to use the Picard
iterative procedure
y(n+1)=f(x)+lambda*Int[G(x,t)*F(t,y(n)),{t,a,b}] to generate
the solution y(x). This procedure will converge (wihout the
need for relaxation) under most instances to the correct
solution when n is made large enough. Even the lower order
approximations starting with y0=f(x) yield reasonable results.
To demonstrate this point consider the boundary value problem
y"+x*y=1 subject to y(0)=y(1)=0. We first converted this to a
Fredholm integral equation containing the Triangular Kernel
T(x,t) and then started the iteration with y0=f(x)=-x*(1-x)/2.
The next iteration y1 is found to already give a very accurate
result. One knows this is so by observing that the next
iteration y2 essentially coincides with y1 in the range
0<x<1. Note that the boundary conditions are satisfied
by all y(n)s. In solving a problem numericaly one can use the
Runge-Kutta method using a trial and error approach
involving a guess for the value of the slope y'(0). From our
y2 iteration we know this slope is approximately y'(0)=-
10433/20169=-0.5175. Charles Emile Picard(1856-1941)was a
professor of mathematics at the Sorbonne in Paris and is best
remembered for his method of successive approximations now
referred to as the Picard Method. He worked on functional
analysis and differential equations and was an outstanding
teacher. Click HERE for
a picture of Picard. You can find pictures of lots of other
mathematicians by going
to- http://www-history.mcs.st-andrews.ac.uk/history/Mathematicians/**

**CONSTRUCTION
AND
GRAPHING OF GREEN'S FUNCTIONS: -We have shown in
class that a linear second order differential operator L and an
associated set of homogeneous boundary conditions has a unique
Greens function associated with it. This Greens function G(x,t)
is a two part function which satisfies LG=0 plus the
requirements that G(x,t) is continous at x=t and that
(dG+/dx)-(dG-/dx)=-1/p(t) , where p(t) is the variable
coefficient multiplying the second derivative in the operator L.
The Greens function will be symmetric whenever L is
self-adjoint. To plot such a two part function using programs
such as Matlab or Maple you simply multiply G- by the window
function H(x-a)-H(x-t) and add to it G+ times the window
function H(x-t)-H(x-b). Here H(x-c) is the Heaviside step
function which has value 0 for x<c and value 1 for x>c. In
the accompanying graph I show G(x,t) corresponding to G"+G=0
with y(0)=y(Inf)=0.**

**George Green(1793-1841) was the son of a baker in Sneinton,
Nottingham, England . He was a mathematical genius, although his
formal education had stopped with grade shool. His occupation
was that of miller and his hobby was mathematics, which he
largely learned on his own. In 1828 he wrote the first of his
ten or so remarkable scientific papers. This first paper,
entitled "An Essay on the Application of Mathematical Analysis
to Theories of Electricity and Magnetism" , brought him to the
attention of the scientific establishment and he was invited and
then began to attend Cambridge University as an undergraduate at
the ripe old age of forty. Unfortunately his health deteriorated
and he died a few years later at the age of 48. His name is
associated with Green's theorem and the various Green's formulas
and he is also credited with invention of the Green's function.
He was never married but had seven children with the daughter of
his mill foreman. I couldn't find any photos of Green.**

**GREENS
FUNCTION FOR A THIRD ORDER OPERATOR: -When converting
boundary value problems governed by an nth order linear ODE one
follows the recipe that L(G)=0 and that all derivatives from
zeroth to n-2 are continous at x=t. Also there is the jump
condition that the (n-1) derivative of G+ minus the (n-1)
derivative of G- equals to -1/p(x) at x=t. Applying this
procedure to G"'=0 with y(0)=y'(0)=y'(1)=0 yields the two part
Greens function shown on the graph. Note that here one does not
see a kink in the function at x=t since the first derivative is
continous. On the other hand the curvature changes.**

**SYMMETRIC KERNELS , THE ORTHOGONALITY OF
Y(X) , AND THE REALNESS OF THE EIGENVALUES:
In class we proved that a homogeneous Fredholm equation
with a symmetric kernel K(x,t)=K(t,x) has real eigenvalues and
its corresponding eigenfunctions are orthogonal. We
demonstrate these results for the equation**

**y(x)=l
Int[cos(x-t)*y(t), t=0..p/2].
Note that here K(x,t) is not only symmetric but is also of the
PG type . One can see by inspection that it has a solution of
the type y(x)=A*cos(x)+B*sin(x) and that there will be two
real eigenvalues. We find that [( p^2/2--4)/16]*l^2-( p/2)*l+1=0 with real roots l1=0.7779 and l2=3.5039
and the corresponding orthonormal eigenfunctions are
Y1=0.6237*[cos(x)+sin(x)] and Y2=1.3236*[cos(x)-sin(x)] .
These are plotted in the accompanying graph.**

**It is clear from these results that
Int(Y1*Y2, x=0..Pi/2) =0 as expected.**

**ORTHOGONALITY
OF THE EIGENFUNCTIONS OF A HOMOGENEOUS FREDHOLM EQUATION: -We
have
shown that a homogeneous Fredholm equation with a symmetric
kernel is equivalent to a Sturm-Louiville differential system.
From this result we know that all eigenvalues of such equations
are real and that their eigenfunctions satisfy the orthogonality
condition Int[r(x)y(n,x)y(m,x),{x,a,b}]=0 if integer n is
unequal to m. Here r(x) is the weight function and represents
the coefficient of x multiplying the eigenvalue lambda in the
original equation. For the case of the Bessel functions of order
zero, this leads to the condition
Int[xJ(0,mun*x)*J(m,mum*x),{x,0,1}]=0, where mun and mum are the
nth and mth zero of J(0,x). The accompanying graph shows the
functions J(0,mun*x).**

**EXPANDING
THE GAUSSIAN IN TERMS OF LEGENDRE POLYNOMIALS: We have found in class that the solutions to
the Sturm-Liouville differential equation (p(x)y'(x))'+[q(x)+ l*r(x)]y(x)=0 are derivable from a
solution of the equivalent homogeneous Fredholm equation with
a symmetric kernel. Thus we know that the eigenfunctions y _{n}(x)
form an orthogonal set with a weight function of r(x) in the
range a<x<b. That is Int[r(x)*y_{n} (x)*y_{m}(x),
x=a..b]=Const.*d_{ m,n}
,where d_{n,m} is the
Kroenecker delta equal to one when n=m and zero when this is
not so. This orthogonality result allows one to take the
functions y_{n}(x and expand any function f(x) in
terms of these as f(x)=Sum[C_{n}*y_{n} (x),
n=0..infinity] in the range a<x<b. The most common set
of functions used are the sin(nx) and cos(nx) functions,
although one can use any other orthogonal set representing
solutions to the S.L. system. We show you here the expansion
of the Gaussian using the Legendre polynomials which are
orthogonal in -1<x<+1 and have r(x)=1 as a weight
function. Note we get a very good approximation using just
three terms of even Legendre polynomials.**

**HILBERT-SCHMIDT SOLUTION OF A
NON-HOMOGENEOUS FREDHOLM EQUATION WITH SYMMETRIC KERNEL:-We
have shown that a homogeneous Fredholm equation with a symmetric
kernel has real eigenvalues and orthogonal eigenfunctions.
We show here how Hilbert and Schmidt were able to use these
properties to solve non-homogeneous Fredholm equations with
symmetric kernels.**

**SOLUTION OF
A NON-LINEAR INTEGRAL EQUATION: -In many instances one can solve homogeneous
non-linear integral equations of the type
y(x)=lambda*int[K(x,t)*F(t,y(t)),{t,a,b}], where F is a
non-linear function of t and y(t), by inspection noting that
the x funtional behaviour must be the same on both sides of
the equation. We demonstrate this for
y(x)=lambda*Int[sin(x-t)*y(t)^2,{t,0,pi/2}]. Here it is clear
that y(x) must be of the form y(x)=A*sin(x)+B*cos(x), with A
and B being constant obtained by solving two simultaneous
algebraic equations. The only non-trivial real roots in this
case occur for A=3 and B=-3 and the solution thus assumes
the periodic form shown in the graph. Note that for such
nonlinear problems the eigenvalue lambda looses its
significance and can take on any value one wishes to set. We
have set it to one. The above equation is also very close in
form to the standard Hammerstein equation , except that the
latter requires a symmetric kernel, which is not the case in
our example.**

**MULTIPLE
SOLUTIONS FOR A NON-LINEAR INTEGRAL EQUATION: An interesting property of non-linear integral
equations is that they can have several solutions, one
solution, or no solutions for the same equation by simply
changing the value of a constant in the equation. We
demonstrate this fact here for the non-linear IE given by
y(x)= l*Int[(x+t)*exp(y(t)),
{t=0..1}], with l being a variable parameter. For
this equation all solutions must have the form y(x)=Ax+B
with the values of A and B determined by simmultaneously
solving the algebraic equations
A=Lambda*Int[exp(At+B),{t=0..1}] and
B=Lambda*Int[t*exp(At+B),{t=0..1}]. A little manipulation
shows these equalities are equivalent to
B=-1+A/[1-exp(-A)] and B=ln[A^2/(Lambda*(exp(A)-1)].
Plotting these last two equalities and looking for their
intersections shows that one has two, one or no solutions
for Lambda <0.332, Lambda=0.332, and Lambda>0.332,
respectively. Similar type of solution behavior occurs for
certain types of non-linear differential equations such as
for the equation of Bratu.**

**DERIVATION OF THE EULER-LAGRANGE EQUATION: -The fundamental problem in the calculus of
variations is to determine the form of a function y(x) which
extremizes an integral involving the functional
F(y(x),y'(x),x). By carrying out a first variation on such an
integral(I=Int[F(y,y',x),{x,A,B}]) one comes to the conclusion
that y(x) must satisfy the Euler-Lagrange equation. We give
this derivation here. Note that a first variation only
guarantees that y(x) produces an extremum of the integral I.
To see whether or not this corresponds to a maximum or a
minimum requires that one look at the sign of the second
variation of I with respect to the parameter epsilon.**

**THE BELTRAMI IDENTITY: The Euler- Lagrange Equation can be integrated
once and hence reduced to a first order ODE for the special
case where the integrant F(x,y,yï¿½) in the variational
integral is a function of y and yï¿½ only. Here , after
multipying the E-L by y', one finds that-**

**y'[(F _{y} )-d/dx
(F_{y'})] = d/dx[F-y'F_{y'}]=0. Thus one concludes that y' F_{ y'} -F=Const. This first order ODE is generally easier
to solve than the correpondng original second order equation
for the case where d/dx=yï¿½F y +yï¿½ F yï¿½. The result
is referred to as the Beltrami Identity(1886) although
Euler had already found it some sixty years earlier. Thus when
F=y+(yï¿½)^2, one can solve either the first order ODE
yï¿½^2-y=C or the second order ODE yï¿½=1/2 to find the same
quadratic curve y(x).**

**BRACHISTOCHRONE
PROBLEM
OF BERNOULLI:-This
cycloidal trajectory gives the minimum time path of a bead
sliding down a curved wire in the absence of friction. First
posed and solved by Johann Bernoulli in 1696. Correct
solutions in a subsequent contest were also given by Leibnitz,
Newton , L'Hospital and Johann's older brother Jacob. The
Bernoullis were a remarkable Swiss family residing in Basel
and producing some eight world renowned mathematicians, the
three best known of which are Jacob Bernoulli(1654-1705),
his younger brother Johann Bernoulli(1667-1748), and
Johann's son Daniel Bernoulli(1700-1782). The fierce disputes
among the Bernoulli brothers , and later between Johann and
his son Daniel, over mathematical priority, are legendary.**

**DERIVATION OF THE
BRACHISTOCHRONE BY VARIATIONAL METHODS: -In the
brachistochrone problem one encounters the first order
non-linear ODE (dy/dx)^2=-(1+y)/y which governs the minimum time
path of a particle sliding down a wire. Go HERE to view a pdf
of the solution which leads to the familiar cycloid curve to
which Bernoulli was referring in his challenge.**

**GOLDSCHMIDT PROBLEM: -One of the classical problems in the calculus
of variations deals with the minimum surface area a soap film
will assume when held between two rings located at +xo and -xo
along the x axis. From the Euler-Lagrange equation one has
here that such a minimum axisymmetric surface satisfies the
equation y'(x)=sqrt[(y/c)^2-1], where x is the axial and y the
radial distance, respectively. A solution of this equation is
the catenoid surface y(x)=c*cosh(x/c) and the ring boundary
conditions requires that p=cosh(p*xo) , with c=1/p. It was
Goldschidt who first observed that this last transcendental
expression can have two solutions, one solution and no
solutions. Of particular interest is the one solution case
which occurs just before the film breaks as xo becomes too
large. This breakpoint is found where the transcendental
equation is satisfied and also the derivative condition
1=xo*sinh(xo*p) is met. A calculation shows this point to
occur at p=1/c=1.81017 and xo=0.66274 A plot of this surface
is shown by going HERE .**

**SOLUTION FOR EXTREMUM SURFACE HAVING
NON-SYMMETRIC END CONDITIONS: The
above discussed Goldschmidt problem dealt with symmetric end
conditions and showed that above certain values for xo no
solution can exist. This fact continues to hold when dealing
with non-symmetric end conditions. We show here a solution for
a surface where y(0)=0.2 and y(0.25)=1. Note that the distance
between the ends is kept at the small value of 0.25 in
order for a solution to be possible. The break point will
occur when the distance x1-x0 between the ends goes to 0.457.**

**MINIMUM TIME PATH BETWEEN TWO POINTS:-One of the
interesting problems encountered in the calculus of variations
is the path a particle or light ray will take in a velocity
dependent medium. The transit time is given by the integral t=Int{ds/v(x,y),{x,A,B]}, with
ds^2=dx^2+dy^2, v(x,y) being the coordinate dependent local
velocity, and A and B values of x at the endpoints of
the path. We show here a set of solutions for the special case
of v=a+b*y in which case the Euler-Lagrange equation for the
functional F(y,y')=sqrt(1+y'^2)/(a+b*y) reduces to the first
order equation 1+y'^2=[-1/(c*(a+b*y)]^2. This equation has the
solution [c*b*x-k]^2+[c*(a+b*y)]^2=1, where c and k are
constants. It will be recognized that this solution represents
a family of circles of radius 1/c*b centered at
[k/(c*b),-a/b]. Note that the extremum path in time tends to
lie predominantly in the high velocity regions. Light rays in
optical fibers follow such curved paths when propagating
through gradient index glass.**

**LIGHT
PATH THROUGH A GLASS SANDWICH:- Consider two equally thick glass plates of
different but constant index of refraction n forming a
sandwich. We assume that n=n _{ 1} for 0<x<0.5
and n=n_{2} for 0.5<x<1. If a light ray
starts at A located at (0,0) and goes to B located at (1,1),
it will take a minimum time path consisting of two straight
lines of different slope which intersect at an interfacial
intersection point (x, 0.5). We show here this point versus
the index of refraction ratio y=n_{ 2} /n_{ 1}.
This result is easiest to derive using Snell's Law of
refraction which (as shown in class) follows from the Fermat
Principle.**

**LIGHT
PATHS IN GRADIENT INDEX GLASS: -We have shown that as a consequence of the
Fermat Principle, light rays in a gradient index glass
generally do not travel in straight lines but rather choose
those paths which minimize the transit time between two end
points A and B. For the special case of gradient index fibers
where the index of refraction n=n(y) is a function of the
transverse direction only, one must solve the equation
1+y'^2=[n(y)/c]^2. We show one such solution for the case
n=sqrt(2)*(1+e*y) in the accompanying graph. Note the big
difference a very small change in the index of refraction has.
For optical fibers one wants to have the index of refraction
decrease with increasing y (ie e<0) so that the light ray
bends back down toward the fiber axis. We point out that these
paths will have no real physical meaning in those regions
where n becomes less than one or where n is much above 1.6.**

**LIGHT RAY PROPAGATION IN A SELFLOC OPTICAL
FIBER: -The Euler-Lagrange
equation for light propagation through an opticakl fiber with
transverse dependent index of refraction n=n(y) can often be
solved in analytic form. One such index of refraction
distribution is the selfloc form n(y)=sqrt[c^2*(1-a^2*y^2)],
where c is the index on the axis and 'a' is a constant which
indiccates the change in n with increasing y. For this
distribution one finds the closed form solution
y=sqrt(b^2-1)*sin(a*b*x)/(a*b), with 1+[dy(0)/dx]^2=b^2. For
the three curves given here the launch angle is 45
degrees but the value of the constant 'a' differs. Note that a
smaller transverse gradient in n(y) lead to larger wavelength
and amplitude for the sinusoidal path.**

GEODESIC ON A CONE AND SPHERE: -For the case of a conical surface the Euler-Lagrange equation leads to the corksrew like geodesic as shown. For a sphere the geodesic is part of the great circle route which appears curved on a Mercator projection but actually represent the shortest distance given by the intesection of a plane passing through the sphere center and the sphere surface containing the end points A and B. Geodesics play a major role in Einstein's General Theory of Relativity.

SHORTEST DISTANCE BETWEEN TWO POINTS ON A SPHERE: We have shown in class that the shortest distance between two points on a sphere is formed by the intersection of the sphere's surface with a plane passing through the sphere center and containing the two points. This is the geodesic for a sphere. Now the question which arises is -what is the distance between two points when one moves along the geodesic? In terms of spherical coordinates one has DIST=R*Int[sqrt(sin(theta)^2 dphi^2+dtheta^2] where R is the earth radius. Since the equation for the geodesic on a sphere was found to be a rather complicated expression, it is impractical to attempt this integration. One rather uses an approach involving spherical trigonometry and a terrestrial triangle having the north pole N, the location of position P1 and location of P2 as its three vertices. By clicking on the title, you see how easy this approach is. We present a specific calculation showing that the shortest distance from London to New York is about 3456 miles.

**NEWTON'S PROJECTILE DRAG PROBLEM: -One
of the better known problems in variational calculus is to
determine the shape of an axisymmetric projectile which
minimizes the air drag. Consider such a body y=y(x) and look at
the axial component of the drag force produced by a single
molecule hitting the surface. For a rarified atmosphere this
force is simply the momentum change relative to the surface
normal of 2mvcos(phi) multiplied by a constant times cos(phi).
Also from the geometry one has that sin(phi)^2=y'/(1+y'^2). This
means the total drag force F=2a*Int[y*y'^3/(1+y'^2),{x,b,c}],
where a is a constant and b and c are axial values of the
projectile end points. The Euler-Lagrange equation shows this
last integral is extremized when y=a*(1+y'^2)^2/(y')^3. We give
the Runge-Kutta numerical solution of this non-linear
differential equation in the accompanying graph for the case
a=1/40, y(0)=0.1 and y'(0)=1. One needs to point out that this
shape( which does look a lot like early rifle bullets) does not
represent the true minimum drag shape when considering drag
produced at normal air pressures. In the latter instance the
shape would need to be time dependent changing with the bullet's
varying Mach and Reynolds numbers during its flight.**

**EFFECTIVE
POTENTIAL OF A PARTICLE SLIDING DOWN A ROTATING HOOP: -We
have
developed
the Lagrange equations for the motion of a particle under
conservative conditions starting with the Principle of Least
Action. One such problem concerns the motion of a particle of
mass m along a hoop of radius 'a' rotating about a
vertical axis with constant angular velocity omega. Here one
finds that the Lagrange equation leads to the conservation of
total energy statement E=(m/2)*[r*d q/dt]^2+V( q), where the effective potential energy
V( q ) equals mga[cos( q)-B*sin^2( q)]
and B=a*( w)^2/(2*g) is the ratio of
the centrifugal to the gravitational acceleration. A plot of V( q) is given in the graph. Notice that at
larger B's one can have a stable point at a point other than the
bottom of the hoop.**

**THE 2D
HARMONIC OSCILLATOR AND THE LISSAJOUS FIGURES: Another problem easily treated by the
Lagrangian approach is that of a mass m suspended by four
springs of equal spring constant k acting along the positive
and negative x and y axis. The Lagrangian here is
L=(m/2)(x'^2+y'^2)-k*(x^2+y^2). The two corresponding Lagrange
equations thus are de-coupled and become (after normalization)
X"+X=0 and Y"+Y=0 with the solutions X=sin(tau) and
Y=cos(tau+C). Here X=x/xmax and Y=y/ymax, and
tau=t*sqrt(2*k/m)+Const. We give the plots of the mass
position in 2D as a function of time for several different
values of the phase C. These figures are the familiar ones of
Lissajous and are also the ones seen with an oscilloscope when
combining two out of phase periodic voltage signals when fed
in along the x and y coordinates.**

**EXTENDED EULER-LAGRANGE EQUATION FOR
PROBLEMS WITH A CONSTRAINT: -We can modify the
standard Euler-Lagrange equation to take into account not only
the extremizing an integral of the type I=Int[F(x,y,y')dx] but
doing so when y(x) is subject to the extra integral constraint
K=Int[G(x,y,y')dx]=Const..This modfication is shown here and
will be developed in more detail during class.**

**CATENARY CURVES: -These
are curves in the shape of hyperbolic cosine functions which
represent the minimum potential energy state of a uniform rope
hanging under its own weight in a constant gravitational
field. The governing differential equation, directly
obtainable from the Euler-Lagrange equation with constraint,
is [(y+b)/c)^2=1+(dy/dx)^2, where b and c are constants. This
equation can be solved to yield y=-b+c*cosh(x/c) for solutions
symmetric about x=0. We show here some catenaries of
different length L tied at points x=1, y=0 and x= -1, y=0 .
Note that, although these curves are reminiscent of the
brachistochrone, they differ slightly in shape from it(click HERE) and the transit time will always be longer
along a catenary than along an equal length
brachistochrone(cycloid).**

**PROBLEM OF
DIDO: -This is the problem
of finding the largest area A which can be enclosed in a curve
of fixed length L and falls under the catagory of
extremizing an integral I subject to constraint K. Here
F=y*sqrt[1+(dy/dx)^2] and G=sqrt[1+(dy/dx)^2] so that the
Euler-Lagrange equation becomes( after one integration)
1+(dy/dx)^2=(lambda)^2/(a-y)^2, where a and lambda are
constants. This equation can be solved in closed form to yield
the off-center circle (x+b)^2+(y-a)^2=(lambda)^2.
Several examples of such circles are given in the accompanying
graph. Some of you will recognize Dido as the queen of
Carthage in Virgil's Aeneid. The story goes that she was
promised as much land by the local ruler in Tunesia as she
could border with a rope made from a single cowhide. She
maximized the land area by laying the rope in a semi-circle
with the Mediterranean sea forming the remaining part of the
border. Hence the problem of Dido, also known as the
Isoperimetric Problem.**

**VARIATIONAL FORM OF THE STURM-LIOUVILLE
EQUATION: -We
can
obtain a variational form for the 2nd order and self-adjoint
Sturm-Liouville differential equation (py')'+[-q+(lambda)r]y=0
subject to a set of homogeneous boundary conditions at x=a and
x=b by multiplying by y and integrating. I show the procedure
here. Once the variational form has been obtained, one can use
it to rapidly calculate an approximate value for the lowest
eigenvalue without knowing the exact form of the corresponding
eigenfunction. I remember using a related variational approch
many years ago in my PhD dissertation in which I needed to
solve certain 6th and 8th order eigenvalue problems in order
to obtain values for the critical Taylor numbers in rotating
MHD flows. Charles Francois Sturm(1803-1855) and Joseph
Liouville(1809-1882) where professors at the Faculte des
Sciences and Ecole Polytechnique, respectively. The latter was
a prolific writer , publishing some 400 scientific papers in
his lifetime while still having time to hold elected office in
the French National Assembly .**

**ESTIMATION
FOR THE LOWEST EIGENVALUE: -We
have shown that one can take the standard
Sturm-Liouville equation [p(x)*y']'+[-q(x)+lambda*r(x)]y=0 and
recast it in the variational form
Int[p(x)*[y'(x)]^2+q(x)*y(x)^2,{x,a,b}]=lambda*Int[r(x)*y(x)^2,{x,a,b}].
Thus one can take any trail function phi(x), satisfying the
homogeneous boundary conditions at the end points x=a and x=b,
to obtain an approximation for the lowest eigenvalue lambda.
We have carried out this procedure for the harmonic oscillator
equation y"+lambda*y=0 subject to y(0)=y(1)==0 using the trail
function phi(x)=[x(1-x)]^n. As is seen in the graph , it
yields as a best approximation the value of lambda=
9.8989 at n=1.112. This result is within 0.3% of the exact
value lambda=Pi^2=9.8696.... and suggests that the form
of phi(x) near this point must be fairly close to the shape of
the exact eigenfunction y(x)=sin(Pi*x). Further improvements
are possible using Galerkin and Rayleigh-Ritz procedures which
also yield information on the higher eigenvalues. Click HERE for
a
second example of an eigenvalue estimation.**

**NUMERICAL
DETERMINATION
OF THE EIGENVALUES OF AN ODE: We have shown in our discussions above
how to obtain estimates for the lowest eigenvalue of boundary
value problems governed by second order differential equations
of the Sturm-Liouville type. To see how accurate these
estimates are we need to compare things with values obtained
by numerical methods. A neat way to find values of lambda
numerically is to first solve a related initial value problem
satisfying the left boundary condition and then stretch the
independent coordinate X until the solution passes through the
right limit of the corresponding boundary value problem. Let
me demonstrate for the equation Y"+lambda*X^2*Y=0, Y(0)=Y(1)=0. Consider first the related initial value
problem y"+x^2*y=0 subject to y(0)=0, y'(0)=1 and solve it
numerically via Runge-Kutta. I show you a graph of this
solution here. Note the zeros occur at x*=[2.358, 3.437,
4.252, 4.736,...]. Now introduce the coordinate streching
x=(x*)X. This leads to the conclusion that
lambda=(x*)^4=[30.92, 137.93, 326.87, 593.61,..]. A
variational estimate using a variational approach for this
problem would thus yield a value for the lowest lambda
slightly larger than 30.92. Note that this procedure works
only as long as r(x)=x^n in the Sturm-Liouville equation. It
would also not work well for higher order boundary value
problems( such as the vibrating beam problem discussed below )
since it requires the guess for several derivatives of the
eigenfunction at the left boundary.**

**LOWEST EIGENVALUE FOR A SIXTH ORDER
EQUATION: -Variational
techniques can be extended to higher order boundary value
problems such as the sixth oder ODE encountered in determining
the lowest rotational speed for the onset of secondary flow in
a viscous fluid between concentric rotating cylinders(Taylor
Problem) or the lowest adverse temperature gradient at the
onset of convection in a horizontal fluid layer(Benard
Problem). Here one has the equation [D^2-a^2]^3*y+T*a^2*y=0
with D=d/dx. The eigenvalue of the problem is T and it varies
with the value of the parameter a. The boundary conditions are
that y, Dy and [D^2-a^2]^2y all vanish at x=0.5 and x=-0.5.
Clearly to find a trial function phi which can satisfy all
these six boundary conditions would be a formadible task.
However, one can get around this difficulty by carrying out a
variational procedure on the two simultaneous equations which
are equivalent to this sixth order form. Such a new procedure
was applied the first time to this problem in my Princeton Phd
dissertation and leads to the variational estimate
T<Int[(u")^2+2a^2*(u')^2
+a^4*u^2,{x,-0.5,0.5}]*Int[(v')^2+a^2v^2]/[a*Int[u*v,{x,-0.5,0.5}]^2,
where u and v are trial functions satisfying the much simpler
boundary conditions u=Du=v=0 at x=0.5 and -0.5. I show you
here an estimate for the eigenvalue T versus the wavenumber
'a' using such a variational approach. The trail functions
u=[1-(2x)^2]^2 and v=1-(2x)^2 have been used. A comparison
with an exact numerical calculation shows that these
variational estimates are very good and will always lead to an
upper bound for T. More details on the convergence of this
variational method can be found HERE.**

**SOLUTION OF THE CLAMPED BEAM EQUATION: -Another
classical equation which can be nicely treated by variational
methods is the clamped beam equation. It reads d ^{4}y/dx^{4}
-ly=0 subjected to the four boundary
conditions y=y'=0 at x=-1 and x=+1. Its variational form is
obtained by multiplying the equation by y, integrating over
-1<x<1, and then integrating by parts. The resultant
variational form reads l = <(y")^{
2}>/<y^{2}>, with the bracket notation
representing the integral over the indicated range of x. Using
the very simple trial function j=(1-x^{2})^{
2} to approximate the lowst eigenfunction y1, one finds
that the lowest eigenvalue must be l<31.5.
A not bad result considering that the exact lowest eigenvalue is
l =31.2852. Improvements in this
upper bound can be found by use of the Rayleigh Ritz method as
we will demonstrate in the upcoming lecture . Note that the
above differential equation lends itself to exact solutions .
These are readily shown to have the even form y_{even}=cosh(kx)/cosh(k)-cos(kx)/cos(k)
and
the odd form y_{odd}=sinh(kx)/sinh(k)-sin(kx)/sin(k),
with the eigenvalues l=k^{4}
determined from the transcendental relations tanh(k)+tan(k)=0
for the even modes and coth(k)-cot(k)=0 for the odd modes. These
beam functions are of considerable utility in certain Galerkin
procedures as a convenient set of orthogonal trial functions
requiring the indicated four boundary conditions(see, for
example, one of my early papers on the "Convective Instability
of a Hydromagnetic Fluid within a Rectangular Cavity", Int. J.
Heat and Mass Transfer 8, 35-41, 1965.)**

**RAYLEIGH-RITZ-GALERKIN
METHOD
FOR OBTAINING IMPROVED EIGENVALUE AND EIGENFUNCTION
ESTIMATES:Consider again the Sturm-Louiville
equation and this time approximate its solution by the
function y(x)=c _{1} j_{1}+
c_{2 }j_{2}+c_{3}j_{3} +--. Here j_{n }represents trial functions
each satisfying the imposed homogeneous boundary conditions of
the problem. If one now multiplies the Sturm- Liouville equation
by j_{m} , substitutes the
above series for y(x) , and then integrates over the range
a<x<b, there results the equality Sum[c_{ n}
(A_{nm}-lB_{ nm}
),{n=1..N}]=0, where A_{nm}=<p(x) j_{ m}(x)' j_{n}(x)'+q(x)j_{n}(x) j_{m}(x)>
and B_{nm}=<r(x) j_{n}(x) j_{m}>. Here < >
indicates integration of the indicated known functions over
a<x<b. If this procedure is repeated with all other N-1
trial functions, there results a set of N linear homogeneous
algebraic equations for the coefficients c_{n .} For
non-trivial solutions for c_{n}, the determinant of the
coefficient matrix of these equations must be zero.
Evaluating this determinant yields N roots for the l s which form upper limits for the first
N eigenvalues of the problem . Once these are found, one can
then evaluate the ratio of the various c_{n}s by
standard algebraic methods involving the minors of the
coefficvients.**

**CHARACTER OF THE EXTREMA OF A TWO-VARIABLE
FUNCTION: -Connected
with
Lagrange multipliers for a problem with no-constraints is the
determination of the character of the extrema of the function
z=f(x,y). Looking at this function as a 2D contour plot within
the x-y plane , one finds that it has an extremum wherever
dz/ds=grad(f(x,y)-z).[icos(theta)+jsin(theta)]=(df/dx)*[cos(theta)]+(df/dy)*[sin(theta)]=0
.
Here the derivatives with respect to x and y represent partial
derivatives. Taking a second derivative one finds at the
extremum that
d^2z/ds^2=[]=(df/dy)^2A-2*(df/dy)*(df/dx)*B+df/dx)^2*C, where
A=d^2f/dx^2, B=d^2f/(dxdy), and C=d^2f/dy^2. This expression
can have no zeros if B^2<AC and thus makes its extrema
there either a maximum or a minimum. A maximum
corresponds to []<0 near the extremum point while a minimum
occurs when []>0. The sign of A and C determine which type
one has. When B^2>AC it is possible for [] to have a zero
at some angle theta so that the critical point changes its
sign as theta is varied. This corresponds to a saddlepoint. We
demonstrate these results by looking at the contour plot of
the function z(x,y)=0.5*x^2-0.25*x^4-0.5*y^2. This
function has three extrema, two being maxima and one a saddle
point.**

**
**SOME ADDITIONAL
MATHEMATICAL TOPICS OF POSSIBLE INTEREST

TOPICS FOLLOW CHRONOLOGICAL ORDER WITH THE LATEST
CONTRIBUTIONS AT THE BOTTOM OF THE PAGE

**DIVERSION ON FRACTALS: Most of you are familiar with the topic of
fractals and how they are generated by certain iterative
relations. There are available on the WEB quite a few programs
which will generate fractals within the complex Z=X+i*Y plane.
A very nice MATLAB program is given by Alberto Strumia of the
University of Bologna and found at http://eulero.ing.unibo.it/. We show you HEREan interesting fractal pattern generated using
a modified version of his program for the special iteration
Z[n+1]=sin(Z[n])-0.085. It exhibits a pleasing spatial and
color symmetry.**

** SHORT DISCOURSE ON WAVELETS: Another very active area of applied
mathematics research is that of wavelets. The technique has
the advantage over windowed fast Fourier transforms(FFTs) by
being able to decompose time varying input data into
simultaneous frequency and time domains . Wavelet approaches
are useful for the filtering of noise from input data and also
for the compression of data for storage.**

**Let me demonstrate a simple application
of wavelets. Consider the input string of data
[F(1),F(2),F(3)....F(N)] , where the Fs represent values of a
time dependent function sampled at N=2^k equally spaced
sub-time intervals n=1, 2...N. This signal can be decomposed
by averaging each of the two neighboring inputs producing a
string of N/2 elements. We call this the first averaging A1.
At the same time, we also obtain a new function(the first
wavelet W1) which measures the difference between F(n) and the
average for the two time unit sub-intervals. Repeating the
procedure by taking the average of the previous average now
over four time units, leads to the second average A2
represented by a string of length N/4. Also taking the
difference of it with the earlier averaged value, one finds a
second wavelet W2. Repeating this procedure until there is
just a single overall average Ak left, one arrives at the
decomposition F=Ak+W1+W2+...+Wk of the original signal. One can now apply all
kinds of modifications to the W strings , such as setting
terms in the Ws which are small to zero. This leads to sparse
matrixes which can be very efficiently stored and manipulated
for later reconstruction. By clicking on WAVELETS you can see the graph of a wavelet
decomposition for the eight element string F=[2, 4, 8, 6, 3,
1, -2, -4], where A1=[3, 7, 2, -3 ], A2=[5, -0.5],
A3=[2.25],W1=[-1, 1, 1, -1, 1, -1, 1, -1], W2=[-2, 2, 2.5,
-2.5], and W3=[2.75, -2.75]. Thus, for example,
F[4]=2.25-1+2+2.75=6. Note the process involves just
addition and subtraction and no integration and is thus well
suited for computer manipulations and involves only O(N)
operations as opposed to O(N ln(N)) for FFTs.**

**SOME MANIPULATIONS INVOLVING THE
GEOMETRIC SERIES: In several of
the lectures above we have made use of the geometric series Sum[x^n, n=0..infinity]=1/(1-x) in deriving things like the Poisson integral
for the circle and for obtaining of some Laurent expansions. I
have recently spent a few hours on manipulating this series to
obtain some other interesting results involving the equality
between certain related infinite series and definite
integrals. You can find some of our results by clicking on the
pdf file GEOMETRIC-SERIES You may want to use these equalities
for obtaining what must be an infinite number of other
identities involving infinite series.**

**ZEROS OF THE RIEMANN ZETA
FUNCTION: One of
the most interesting functions encountered in classical
analysis is the Riemann Zeta funtion defined as
Zeta(x+I*y)=Sum[1/n^(x+I*y), n=1..infinity]. This function has
values Zeta(1+I*0)=infinity, Zeta(2+I*0)=Pi^2/6 and, according
to Riemann's Hypothesis(1863), has zeros only when the real
part x has a value of 1/2. Also, as shown quite early by
Euler, using a similar argument to that used in summing the
geometric series, one has the equality Zeta(z=x+I*y)= 1/
Infinite Product[1-(1/p^z)], where p are the prime
numbers 2, 3, 5, 7, 11,etc.. A recent book "The Prime
Obsession"by John Derbyshire gives an excellent discussion of
attempts over the past 140 years to prove that x must
necessarily be 1/2 for a zero of the Zeta function to exist.
So far no successful proof has been found and it remains one
of the unsolved of the 23 problems posed by Hilbert at the
Second International Congress of Mathematicians at Paris in
1900. That the hypothesis is highly likely to be correct
can be shown by a simple computer calculation as I have
carried out here. Note that when x departs from x=1/2 the
curve for abs[zeta(x+I*y)] no longer reaches zero. An
alternative way to determine the values of y for which the
x=0.5 Zeta function has zeros is to plot the functions
Re[Zeta(0.5+I*y)] and Im[Zeta(0.5+I*y) and then see where both
vanish simultaneously. You can see this procedure by going
to ZETA-ZEROS,
with the blue circles indicating the zero locations.**

**THE EULER-MASCHERONI CONSTANT AND THE
HARMONIC SERIES: In the
derivation of theBessel Function of the second kind one
encounters the constant
gamma=0.577215664901532860606512...comonly referred to as the
Euler-Mascheroni constant. It is defined as the difference
between the divergent harmonic series and the logarihm of
infinity. Go to EULER-MASCHERONI-CONSTANTto to see a more detailed discussion on
this constant and its the relation to the harmonic series with
a finite number of terms and to the digamma
function. In case you are wondering how the mathematics
professor Lorenzo Macheroni(1750-1800) of the University of
Pavia got his name attached to Euler's constant (1735), it
comes from the fact that in 1790 Mascheroni calculated gamma
to 32 places(the last thirteen of which were wrong) while
Euler had only expanded things out to the first sixteen
places.**

**FORMULAS FOR THE SUMMMATION OF INTEGERS
TAKEN TO THE PTH POWER: In
numerous analytical manipulations including those encountered
in calculus and in statistics, it is of interest to sum the
first N integers taken to differrent powers p. The simplest of
such summations is 1+2+3+....+N which most readers will
recognize as equal to N(N+1)/2. What is perhaps not as obvious
is that the sum of the odd intergers equals
1+3+5+...+N=(N+1)^2 and that
1+2^3+3^3+4^3+..+N^3=(1/4)[N(N+1)]^2. Click on PTH-POWER to
see how formulas for the sum 1+2^p+3^p+..+N^p for p=1,
2, 3, 4 and 5 are obtained.**

**IDENTITIES INVOLVING N! : In many of the calculations above we have
encountered the factorial function n! and variations thereof.
Such variations occur in expressing the hypergeometric series
and obtaining the binomial coefficient and also in many other
applications. By going to FACTORIAL-ZEROS you will find a few of the more common
factorial identities including some not found in handbooks.**

**BINOMIAL
THEOREM, PASCAL TRIANGLE, AND THE GAUSSIAN: It is well
known that (a+b)^n can be expanded as an (n+1) order
polynomial involving powers of a and b. Newton was the first
to present a mathematical formula for such a binomial
expansion although the binomial coefficients in the expansion
were known earlier in the form of the Pascal triangle. Also
one notes that the magnitude of the coefficients in a
binomial expansion have the form of a Gaussian distribution
as n approaches infinity. Go to BINOMIAL-THEOREM to see some discussion of these facts and how
they are related to each other.**

**BINARY
NUMBER
SYSTEM: As you are all
aware, we have, during the course of the last thirty years,
entered the so-called digital age in which most of our
information transmission , manipulation and storage is
accomplished in binary form. Be it by means of DVDs, compact
discs, optical cable, computers, the internet, or cell phones,
information is being handled strictly in units of bits, the
basic binary unit. The mathematics underlying this binary
technology dates back to G.W.Leibnitz in the late sixteenth
century and is based on the binary
number system in which there are
just two basic states of on (1) or off (0) . Such binary
states are ideally suited for electronic or optical circuits
and do not suffer the degradation of analogue devices of the
past. By means of such a binary system one can encode,
store, manipulate, and transmit anything , be it
numbers, words, pictures, or music. I present at BINARY a brief discussion of the binary number
system.**

**SUMMATION
OF SERIES USING LAPLACE TRANSFORMS: It is possible to sum many infinite series by
a technique involving the Laplace transform of functions. The
method allows one to sum infinite series from n=0 or n=1 to
n=infinity, by evaluating certain indefinite integrals of the
form Int[ f(t)/(1-exp(-t)), t=0..infinity]. Go to LAPLACE-TRANSFORMS to see some examples.**

**CATALAN
CONSTANT:This constant
named after the Belgium mathematician Eugene
Catalan(1814-1894) has the form G=1-1/9+1/25-1/49+.. The
series converges to the (as yet unproved) irrational number
0.915965....There are numerous other series and integral
representations for G. I develop some of thes ate CATALAN.**

** SERIES BY COMPLEX VARIABLE
METHODS: A final method for
obtaining the sum of infinite series is by a complex variable
approach involving the functions f(z)cot(pi z) and f(z)csc(pi
z) to sum infinite series of the form Sum[f(n),
n=-infinity to +infinity]. We give you a brief summary of the
technique at SERIES-BY-COMPLEX-METHODS-.
**

**POLYGONAL WEB NUMBERS: It is well known that there
are an infinite possible sequences of numbers which can be
generated from the formula F(n)=an^2+b^n+c where a, b
, and c are specified constants. Perhaps one of the best
known of these is the sequence T(n)=[1, 3, 6, 10, 15,
20,...] which is generated by the formula T(n)=n(n+1)/2 and
represents the sequence of numbers representing the
sum of the first n integers. We consider here a variation on
this classical sequence by asking for the number of points
out to the nth row of a polygonal web as shown at POLYGONAL-NUMBERS
Such webs are constructed by drawing radial lines out
from the center of a regular polygon and then constructing
rows at distances out equal to an integer multiple of the
distance from the polygon center to a vertex. Placing
equally spaced points on the various rows yields a
point number represented by the formula F(s,n)=s[n(n-1)+2]/2, where s is
the number of sides of the regular polygon web and n the row
number one has moved out to. This formula is constructed by
recognizing that each sector of the polygon net just
contains the triangular number of points T(n) and thus
the total number of points must be sT(n) minus the number
overlaping points which equals s(n-1)-(s-1). The hexagonal
web numbers(shown as red points in the figure) are
F(6,n)=3n^2-3n+1= [1,7,19,37,61, 91, 127,..]. We can also
consider summing up the reciprocal of the polygonal web
numbers as an infinite series. A little
manipulation(using the Laplace transform method discussed
above, allows one to arrive at the identity-**

**
S(s)=Sum[1/F(s,n),n=1..inf.]=(2/s)*Sum[1/((n-1/2)^2+(2/s-1/4)),n=1..inf.]=
[
4/sqrt(s(8-s))]
*Int[sin[sqrt((8-s)/s)*x/sinh(x),x=0..inf.]
**

**Using this formula for
the case of an octagonal web number sequence , one finds
that
**

**
S(8)=Sum[1/(4n^2-4n+1),n=1..inf.]=0.5*Int[x/sinh(x),x=0..inf.]=(pi^2)/8=1.233700...
**

**GENERATION OF ARCTAN
FORMULAS FOR PI: For
nearly 300 years the standard approach to finding
evermore accurate values for Pi has been the use of
arctan formulas. Although this approach has been
superseded in recent years by the use of
AGM(algebraic-geometric mean) methods applied to the
numerical evaluation of complete elliptic integrals
using the latest supercomputers, it is nevertheless
worthwhile, in view of the countless hours spent by
mathematicians on finding improved version of arctan
formulas, to see how arctan formulas are derived , to
look at those of historical significance, and those
which converge rapidly. Among the more rapidly
converging series is our own four-term arctan formula
for Pi-
**

**
Click on ARXTAN-FORMULA
to see the discussion .
**

**POLYNOMIAL QUOTIENT
APPROXIMATIONS FOR ARCTAN AND ARCSIN: I have recently been experimenting
with obtaining some new approximations for arctan(1/N)
and arcsin(1/N) when N>>1. Using an infinite
integral representation for these functions and then
applying the standard Leibnitz technique one comes up
with some interesting polynomial quotients in N which
lead to excellent approximations for these functions
without ever needing to explicitly evaluate any
integrals. Click on
POLY-QUOTIENTSto
see the development.
**

**GOLDEN
RATIO AND FIBONACCI SEQUENCES: Solutions
of the quadratic equation x(x-1)=1 and the
difference equation F[n+1]=F[n]+F[n-1] are known to
lead to the golden ratio x= (1+sqrt(5))/2=
1.61803... which also equals the ratio of F[n+1]/F[n]
as n goes to infinity. We present at
GOLDEN-RATIO-FIBONACCI some discussion
of these properties and also introduce some additional
Fibonacci like sequences.**

**GAMMA, BETA, AND DIGAMMA
FUNCTIONS: In the above lectures we
encountered numerous functions defined by definite
integrals. These included among others the Bessel
functions, error function, gamma function, and
beta function. We give you at HGDIGAMMA
a more in depth discussion of three of these
functions, namely, the
GAMMA,
BETA, and DIGAMMA function.
**

**FRESNEL INTEGRALS AND THE CORNU
SPIRAL: In our earlier discussion of
solutions to the Helmholtz equation, we considered
the problem of light diffraction by a slit. In
extending such a discussion to diffraction by
circular aperture, one encounters a complex
integral of the type Int[exp-i(ct ^{2}),t=0..x]
which
leads to well the known Fresnel integrals C(x) and
S(x). We discuss at CORNU-FRESNEL
some of their properties.
**

**NUMERICAL EVALUATION OF PI
USING A FOUR TERM ARCTAN FORMULA: We have shown in
earlier discussions above that arctan formulas
for Pi will be rapidly convergent whenever the
values of N appearing in arctan(1/N) are
all large. With this point in mind, we present
at FOUR-TERM-ARCTAN-FOR-Pi
the numerical method we have used to obtain
approximations for Pi . The iteration method
used produces about six additional places of
accuracy in Pi for each additional iteration
and involves only multiplication and division
of integers.
**

**CONTINUED
FRACTIONS: It is possible to express
both rational and irrational numbers plus
functions of x as continued fractions. We
discuss at COINTIUED-FRACTIONS
how this is done, give some of the
better known continued fractions, and show how
they may be used to give good approximations
to certain classical constants such as Pi.
**

**INFINITE
PRODUCTS:
An infinite number of variables and
constants can be expressed as infinite
products and these products can in turn be
used to approximate these quantities to
any desired order of accuracy. We present
at INFINITE-PRODUCTS
some of the better known infinite
products and show how some of these are
derived.
**

**PRIME
NUMBERS: The integers
2, 3, 5, 7, 11, 13,... referred to as the
prime numbers have fascinated
mathematicians for several thousand years.
They are defined as those positive
integers different from one which can be
factored only by themselves and by
one. Several important theorems
concerning these numbers exist and certain
conjectures such as those of Goldbach have
yet to be proven. In the last few decades
the advent of high speed computers have
led to a game of who can find the largest
prime number. In this regard it has been
especially useful to deal with Mersenne
primes as they lend themselves to
relatively easy, but time consuming,
checks. One of the latest Mersenne primes
which has been verified by computer is
the 9.8 million place number (2 ^{32582657}-1)
.
We discuss at PRIMES
some properties of prime numbers and
discuss how one might obtain even larger
primes by a known extention of the
Mersenne postulate that 2^{p}
-1 is prime for certain prime numbers. So
far the first 44 Mersenne primes have been
found. A finacial award of 100K is posted
on the internet for anyone finding the
first p with more than ten million digets.
In playing around with prime numbers, I
noticed that the complex function F(n)= n
exp(i Pi n/4) places all prime
numbers along the intersection of the
Archimedes spiral defined by this F in the
complex plane and the 45 degree diagonal
lines(see
attached PRIME-NUMBER-PATTERN.jpg).
Perhaps such a geometrical interpretation
will in the future be of aid in
determining the prime or non-primeness of
large integers without requiring extensive
divisions by smaller prime numbers.
Finally I give you at GENERATING-PRIMES
a discussion in pdf form of
how to generate large prime and large
composite numbers.**

MORPHING OF THE ULAM SPIRAL: In looking further into ways to graph the location of prime numbers, I recently ran into the the concept of the Ulam Spiral(see http://en.wikipedia.org/wiki/Ulam_spiral ) and some of the elaborate graphings which has been carried out to supposedly show that their computer generated patterns have deeper meaning. To me it seems these patterns simply arise from the fact that primes above n=2 are odd numbers. They say nothing about whether or not a given odd number is prime. Go to MORPHING-ULAM to see our discussion on this topic.

SUMMATION OF INFINITE SERIES USING SPECIAL VALUES OF KNOWN FUNCTIONS: Most analytic functions have infinite series expansions. We show you at SUMMATION how evaluating such series at special values of x yield exact values for certain infinite series. For example, the infinite sum Sum[(-1)^n/(2n+1)!,n=0..infinity]=1-1/3!+1/5!-1/7!+...=j

PROPERTIES OF THE FUNCTION F(a,b)=Sum[1/n!^(a+ib),n=1..infinity): In the previous discussion we examined the values of certain infinite series which are expressible in terms of known functions at fixed values of x. In the process we ran across two infinite series, namely, Sum[1/n!,n=0..infinity] and Sum[1/n!^2,n=1..infinity, for which exact values are given in terms of exp(1)-1 and I(0,2)-1. Here we extend our discussion by looking at the general properties of the new function F(a,b)=Sum[ 1/n!^(a+ib),n=1..infinity] of which the two earlier sums are special cases. Function F(a,b) is reminiscent of the standard Riemann function Zeta(a,b)=Sum[1/n^a+ib,n=1..infinity but differs from it that n inters the sum as factorials. This fact accelerates the convergence rate of the defining series. Go to FUNCTION-VALUE to see some of its properties including the fact that, unlike the Riemann Zeta function, F(a,b) appears to have no zeros for the range of a and b investigated.

LAMBERT FUNCTION AND EQUATION: The first order ODE W'=exp(-W)/(1+W) subject to W(0)=0 has the implicit solution W(x)exp(W(x))=x, where W(x) is known as the Lambert function. Go to LAMBERT-FUNCTION to see some of its properties. Johann Heinrich Lambert(1728-1777) was a polymath born in Alsace who made contributions to mathematics, physics and philosophy. Was a member of the Prussian Academy of Sciences in Berlin where Leonard Euler was his contemporary. He was the first person to prove the irrationality of Pi.

SOLUTION OF INTEGRALS CONTAINING THE GAUSSIAN: In our earlier discussions on integral solutions to the heat conduction and Laplace equations plus the discussion of the saddle point technique, we encountered integrals with infinite limits involving the Gaussian exp-(x^2) in either explicit or implicit form. We derive and evaluate at GAUSSIAN-INTEGRALS some of these integrals in closed form. It is shown , for example, that int[exp(-x^2)/sqrt(x), x=0..infinity]=(1/2)*Gamma(1/4)=1.8128049... . and erfc(b)=(2/sqrt(Pi)*int(exp-(x+b)^2, x=0..infinity).

SOLUTION OF INTEGRALS USING CONTOUR INTEGRATION:Many definite integrals involving infinity as one of its limits can be sloved by contour integration using the Cauchy Theorem and Integral. We present at CONTOUR-INTEGRATION some examples of how this type of integration is carried out. It involves integrating about sinular points and branch points within appropriate chosen closed contours.

BERNOULLI AND EULER POLYNOMIALS AND NUMBERS: There are a set of numbers and functions which arose in the 18th century in connection with the summing of powers of the first N integers. In particular we are talking about the Bernoulli polynomials B

SOLUTION OF AN INFINITE SERIES: We consider the infinite series F(m)=sum[n^m/ 2^n, n=1..infinity] when m is a positive integer and determine the value of F(m) by both direct summation and by evaluation of an equivalent integral. Notice this series will converge for all finite positive values of m since by the ratio test [(n+1)/n]^m<2 as n approaches infinity. The values of F(m) will be found to equal an ever increasing even integers as m increases. Go to INFINITE-SERIESW-SUM to see the discussion.

HEAVISIDE, DIRAC, AND STAIRCASE FUNCTIONS: In numerous applications in analysis we have encountered discontinous functions of which the best known are those of Heaviside and Dirac and functions representing combimation of these. We discuss at HEAVISIDE some of their properties and in particular show how their derivatives are handled by using continous functional representations and then looking at these in the limit as a parameter epsilon approaches zero.

FURTHER PROPERTIES OF COMPLETE ELLIPTIC INTEGRALS: In our earlier discussions on non-linear differential equations we considered in detail the equation (theta)"+sin(theta)=0 representing the swing of a simple pendulum. In solving this equation for time as a function of theta, we encountered the Complete Elliptic Integral K(m). We discuss at ELLIPTIC-INTEGRALS some of the properties of K(m) and its related function E(m) showing that these functions can take a variety of different integral forms. Also we demonstrate that the constant Pi is given exactly by d(K(m)

PROPERTIES OF ARCTAN(Z): In many of the discussions above we encountered the integral int[1/(N^2+z^2), z=a..b) which was shown to have the closed form solution (1/N)[arctan(b/N)-arctan(a/N)]. We discuss at ARCTAN-PROPERTIES some of the properties arctan(z) and show here how this functions arises in several different areas.

EVALUATION OF SOME FINITE SERIES: In most of our above discussions involving the summation of series we have looked at only infinite types such as those defining Airy, Bessel, Hyperbolic functions etc.. We present at FINITE-SERIES a brief summary of how one can sum finite series. Such series are important in determining the error encountered in approximating a function by looking at only the first N terms in their infinite series representations. For example

FINDING LOGARITHMS OF NUMBERS: One can define any number N as b

BINARY REPRESENTATIONS OF PRIME NUMBERS AND CALCULATIONS FOR PRIMALITY: Using the polar diagram for prime numbers developed above, we express the binary form of the primes and use their last three digit patterns to see along which diagonal they fall in our polar diagram. We then show how to evaluate a given number for primality using a modified version of the sieve of Eratosthenes. Go to PRIMES-IN-BINARYfor these discussions.

MODULAR ARITHMETIC FOR THE ADDITION AND MULTIPLICATION OF INTEGERS: Recently we have come up with a way to represent all positive integers graphically using an Archimedes spiral and a set of eight radial lines. The integers are represented as the resultant intersection points with the odd integers lying along diagonal lines and the even integers along the x and y axis. We apply at MODULAR-MATH some of the rules of addition and multiplication they obey.

TRILATERATION AND GPS TECHNOLOGY: Recent years have seen an explosive use of Global Positioning Technology(GPS). It allows one to pinpoint any point on earth at a precise time using signals sent from some 24 to 31 satellites orbiting the earth at an altitude of 12,000 miles and orbital period of 12 hours. To locate a point P(x,y,z,t) will generally require four transmitting satellites within a line of sight and this is achievable with the present total number of satellites for any point on earth. The mathematics behind the x, y, z determination is referred to as trilateration and we give you at TRILATERATION a quick development.

POWERS OF COMPLEX NUMBERS: One knows that a complex number N=a+ib can be conveniently represented as a point in the complex(Argand)plane.Go to COMPLEX-POWERS to see the paths such points will describe as a given complex number is taken to powers n with n extending in a continous manner from n=1 to n=infinity. Among other findings we observe that the paths generally represent exponential spirals and for a special case reduce to unit radius circles. It is shown that i

SHORTEST DISTANCE BETWEEN TWO POINTS ON A SPHERE: In our discussion of the Euler-Lagrange equation in the variational methods section of the above mathematics courses, we encountered the concept of geodesics and evaluated them for specific surfaces such as for planes, cylinders, cones and spheres. We extend these discussions at GEODESICS by looking in more detail at the great circle geodesics on the surface of a sphere. The problem is approached by use of both variational methods and also by spherical geometry. Several examples of distances between points A and B on a sphere are given.

LAGRANGE-MULTIPLIERS: A method to extremize a function F subjected to n constraints is that of Lagrange-Multipliers. We give you at LAGRANGE-MULTIPLIERS a more detailed dicussion of the technique then found above in the Integral Equations and Variational Method Course. Among the examples discussed is the shape of a triangle containing maximum area for a fixed parameter. Also we look at the problem of the shortest distance between two curves.

SPHERICAL HARMONICS: The solution to the 3D Laplace equation when expressed in spherical coordinates reads

V=(Ar

PROPERTIES OF THE COMPOSITE NUMBER N(n,m)=2^(2n+1)+1+6m: In the above disussions on prime numbers we examined the number 2^(2n+1)+1 and found it to always be divisible by three and hence non-prime. Since that time we have come up with a generalization in form of the number N(n,m)=2^(2n+1)+1+6m . We examine at COMPOSITE-NUMBER the properties of this modified number, show it to always be composite, and indicate how it may be used to generate an infinite number of large primes. One such large prime found near the limit of our MAPLE program's capability is the 1303 digit number [2^4327+1+6(282)]/3. See if you can verify that this number is indeed prime.

GOLDEN RATIO: The number Phi=[1+sqrt(5)]/2=1.61803398874989.... represents the Golden Ratio for the ratio of height to width of an ideal rectangular structure according to the early Greeks. Be that as it may, Phi has many interesting mathematical properties. We discuss some of these properrties at GOLDEN-RATIO and show, for example, that Sum[1/Phi

PROPERTIES OF THE EXPONENTIAL FUNCTION F(x)=exp(x): We examine the elementary function exp(x), show how it is derived, give its continous fraction expansion, show its relation to trignometric and hyperbolic functions, discuss the exponential spiral, and use it to rapidly calulate the logarithms of numbers using Newton-Raphson iteration. Go to EXP(X) for the details.

DERIVATION OF THE BAUER FORMULA: An important application of the Helmholtz equation occurs when examining the scattering of an incoming wave by an atomic nucleus. An important identity used in such studies is the Bauer Formula which follows on solving the Helmholtz equation in spherical coordinates and equating the solution to an incoming plane wave Psi=exp(ikz). Go to BAUER-FORMULA to see the details of the derivation.

PROPERTIES OF THE LEGENDRE POLYNOMIALS: An important set of orthogonal polynomials arising in connection with solutions to the Laplace and Helmholtz equations in spherical coordinates are the Legendre polynomials P

SOLUTION OF THE INTEGRAL Int[P

A METHOD FOR QUICKLY FINDING VALUES OF ARCTAN(1/a) FOR LARGE a: In the previous note we examined integrals of the type Int[P

PROPERTIES OF PYRAMIDS: One of the more interesting solids is the standard pyramid of square base and equalateral sides such as exemplified by the pyramids at Giza, Egypt. Although most of you have been through calculating some of the properties of such pyramids starting in your high school calculus course and extending through Lagrange multiplier as done in the above variational methods course, it is of interest to summarize some of the more important properties of pyramids and show how these are derived. Go to PYRAMIDS to see the discussion.

3D INTEGRATION METHODS: Elementary integration as one learns in calculus is predominently confined to one and two dimensional situations with relatively little time spent on full three dimensional calculations.We give you at 3D-INTEGRATION some of the techniques used to handle 3D integrations. In many of the instances we also display 3D computer plots of the volumes one is integrating. Cartesian, Cylindrical, and Spherical Coordinates are used depending on the needs of bounding surfaces encountered.

PROPERTIES OF DUERERS'S CUBE: In the engraving shown in the leadin of this WEB page one sees a figure containing a peculiar looking cube which has two tetrahedrons sliced off of opposing corners. This solid is known as Duerer's Solid or , as I refer to it, as Duerer's Cube. We look at some of the properties of this solid at DUERER-CUBE including 3D graphs showing both the Duerer Cube and its modification.

THE ICOSAHEDRON: One of the more interesting and most complicated of the five Plato polyhedra structures is the Icosahedron. It is characterized by a surface composed of twenty equal area equilateral triangles whose twelve vertexes all lie at a constant distance from its center. This platonic solid has found numerous applications including Geodesic Domes, in the description of Bucky Balls, and even is associated with the structure of certain viruses. Go to ICOSAHEDRON to see our dicussion on this topic.

RADIUS OF CURVATURE, EVOLUTES, AND PROPERTIES OF ELLIPSES: One of the more difficult concepts encountered by students in elementary differential calculus is the radius of curvature of a curve y=f(x) and the generation of its evolute. We give you at EVOLUTES a brief discussion of these comcepts plus discuss some extentions including the focusing capabilities of ellipses.

THE DODECAHEDRON: Another of the five platonic regular solids is the dodecahedron. This solid has twelve regular pentagon faces(F=12), thirty edges(E=30) and 20 vertexes(V=20) and obeys the Euler formula V+F-E=2. We derive at DODECAHEDRON its surface area and volume plus indicate how one may use such a configuration in conjunction with twelve semi-ellipsoids to form an efficient implosion device.

MATHEMATICAL ANALYSIS OF AN AXICON CONCENTRATOR FOR PHOTOVOLTAIC APPLICATIONS: One of the difficuties with converting solar energy into electricity is the high cost of photovoltaic cells. Lately there has been a re-emergence of interest in using solar concentrators for reducing the cost and increasing the efficiency of such conversion processes. We discuss at AXICON the mathematics of an old idea of ours(circa 1980) of using axicon concentrators for uniformly illuminating a conical array of solar cells. The advent of flexible sheets of photovoltaic cells now make this possible although (to our knowledge) no one has actually constructed such a device.

THE GAUSSIAN , BINOMIAL COEFFICIENT, AND THE PASCAL TRIANGLE: It is well known that elements of the Pascal Triangle can be used to expand functions of the form (a+b)^n and that the elements in a given row resemble a Gaussian as the row number n gets large. We give you at PASCAL-GAUSSIAN a discussion of this problem and show that exp-[2k^2/n] can be approximated by (n/2)!^2/{[(n/2+k)!(n/2-k)!] for n/2>k and n->infinity.

2010

THE PROBLEM OF APOLLONIUS: Finding the largest circle which can be placed into the three-sides area formed by three different diameter circles tangent to each other is known as the Problem of Apollonius. Descartes gave an early analytic solution to this problem back in the 1643 and interest in variations including the use of complex numbers have remained in vogue to present day. We discuss this problem at APOLLONIUS and offer some variations.

METHOD FOR QUICKLY EVALUATING ln[1+(1/a)] FOR a>>1: It is possible the find values of ln(x) near x=1 by a technique we have used succesfully earlier for finding arctan(1/a) for large a. The procedure for an accurate numerical determination of the natural logaritm of [1+(1/a)] involves the evaluation of a quotient involving the (2n+1) th order Legendre polynomial for large n and the linear term (x+a). We present the details of such calulations at LOG(x) and show, among other things, that ln([1+(1/a)] = [1/(6a)][(8+15a-30a

KEPLER'S SPHERE STACKING PROBLEM: About 1611 Johannes Kepler, of astronomy fame, conjectured that a certain face-centered form of stacked unit radius spheres left the minimum value of voids possible. This conjecture has just recently been shown to be correct compared to all possible other configurations although earlier mathematicians , such as K. Gauss, already showed it to be correct for lattice based arrangements. We present at SPHERE-STACKING a discussion of the problem and show how one obtains the sphere density value of Pi/(3*sqrt(2).

THE NESTING PROBLEM: Similar to what is encountered with Russian wooden nesting dolls, we look at the problem of continually placing ever smaller figures inside an identically shaped figure of larger size. This is the problem of nesting or embedding and will be looked at NESTING for the special case of circles and spheres. Formulas connecting the dimensions of the smaller figures to the larger are presented. Also explicit values for packing densities are presented.

PRIMES GENERATED BY N[n]=2n(n+1)+1: There are many formulas which can be used to generate prime numbers. One of the best known is the Mersenne Prime formula N[n]=2^n-1. So far only 47 values of n have been found for which N[n] will be prime. We point out that there are many other generating formulas for prime numbers much richer in their prime density. One of these, which we have come up with recently, is the elementary formula N[n]=n^2 +(n+1)^2=2n(n+1)+1 representing the sum of the square of an integer plus the square of its neighbor. Go to PRIMES-NY-FORMULA to see some of its properties including the fact that the high density of primes it produces all end in either 1 or 3.

DERIVATION OF THE CAUSTIC FOR A CIRCLE: One of the more interesting optics problems soluable by ray tracing methods is the determination of the caustic formed by the reflection of light rays from a cylindrical surface. We demonstrate at CAUSTIC how this curve is generated. It is shown that the caustic equals a one cusp epicycloid which, upon an axis transformation, is equivalent to a standard cardioid.

PROPERTIES OF THE NORMAL(OR GAUSSIAN)DISTRIBUTION: We examine the properties in detail at GAUSSIAN and show its approximate relation to the Probability Function. Areas under this curve betwween plus and minus n times the standard deviation are expressed in terms of the error function. The origin of the 68-95-99.7 Rule in statistics is also demonstrated.

INTEGRALS INVOLVING BESSEL FUNCTIONS OF THE FIRST KIND: We have discussed in detail the properties of Bessel Functions earlier. At BESSEL-J(n,x) we look more specifically at methods which can be used to evaluate integrals involving these functions. Laplace Transform techniques, Bessel Function Identities, and Leibnitz Rule applications are found to be useful. A Fourier Series expansion for a parabola using the Bessel Function of the First Kind of order Zero is given.

TETRATION OF THE NUMBER N=a+Ib: We examine the iteration A[n+1]=N^A[n] where A[0]=N=a+Ib is a complex number. This iteration is equivalent to the non-associative power tower N^(N^(N^(N^N))) etc. We show at TETRATION under what conditions the sequencce A[1], A[2], A[3],...converges to a finite value A[infinity] or diverges. A ratio based on the Lambert Function is found to quickly determine A[infinity]. We show that A[n+1]=(sqrt(2))^A[n] with A[0]=sqrt(2) converges to A[infinity]=2 and A[n+1]=I^A[n] with A[0]=I yields A[infinity]=0.4382..+i0.3605.. .

PROPERTIES OF SPIRALS: We examine the various known forms of two-dimentional spirals at DIFFERENT-SPIRALS. After looking at the Archimedes and Logarithmic spirals, we briefly consider other forms such as those of Fermat and Cornu and also indicate how one can create less known spirals. Some attention is also given to spirals which have discontinuities in their derivatives at points along their length.

PADE APPROXIMATES: We show how functions F(x) represented by infiniote series can be approximated by a quotient of two polynomials. The technique, known as the Pade Approximation Method, involves the evaluation of the coefficients of these polynomials such that the first M+L terms of the MacLaurin series for F(x) are satisfied. We give several examples of Pade Approximates, show how one can use canned math programs such as MAPLE to generate them, and indicate how these quotients may be used to determine the location of the singularities of the function F(x). Go to PADE-APPROX for details of these discussions.

TECHNIQUES FOR EVALUATING CERTAIN DEFINITE INTEGRALS: We show at EVAL-INTEGRALS several different methods which may be used to evaluate some non-elementary definite integrals using series expansions, variable subsitutions, Laplace transforms, and integration by parts. Also a method to find approximations for arctan(1/a) and ln(1+1/a) involving integrals containing even Legendre polynomials are presented. We show, for example, that int(sin(ax)*exp(-bx)/x, x=0..infinity) equals arctan(b/a).

GENERATING SOLUTIONS TO THE BIHARMONIC EQUATION BY COMPLEX VARIABLE METHODS: In analogy to the use of complex variable methods to generate a complex velocity potential for 2D inviscid flows, we show here that F(z)=zG(z)+zbarH(z)+J(z)+K(z) generates an infinite number of biharmonic streamfunctions Psi=Re[F(z)] and Psi=Im[F(z)] provided that G(z), H(z), J(z), and K(z) are harmonic functions. Go to BIHARMONIC to see a discussion of the properties of several different F(z)s. The complex function F(z)=zbar/z is found to represent low Reynolds number flow out of a corner between two right angled walls.

MATHEMATICS OF THE BLACK BODY RADIATION LAWS: One of the greatest contribution to classical physics and the beginning of quantum mechanics was the discovery of the black body radiation formula by Max Planck back in 1900. For the first time he showed how one could account for the entire spectrum of radiation emanating from a black body such as the sun. We show you at RADIATION-LAWS the mathematical manipulations involved in formulating Plank's Radiation Law and how it verifies the Wien's Displacement Law and the Stefan-Boltzmann Law.

EUCLIDIAN ALGORITHM AND THE GCD OF NUMBERS: We discuss the Euclidian algorithm for finding the largest common denominator (GCD)between two integers N and M and use the result and the properties of modular arithmetic to find integer solutions of certain Diophantine equations. It is also shown, for example, that integer solutions of x^2=y^3 are found at x=(9/4)(N+1)^2-2 and y=(3/2)(N+1)(x-1)+1 with N=0,ï¿½1,ï¿½2,ï¿½3, etc. The function F(n)=27n^3+81n^2+72n+19 is found to be an excellent generator of prime numbers. Go to EUCLID for the general discussion on the gcd.

PROPRETIES OF ELLIPTIC CURVES AND THEIR USE IN FACTORING LARGE NUMBERS:Curves defined by y^2=x^3+ax+b are known as elliptic curves and historically arose in connection with elliptic integrals. We discuss at PROPERTIES-ELLIPTIC-CURVES some of their characteristics including their ability to aid in the factoring large numbers. Graphs of these elliptic curves for two different values of a and b are presented. The governing equation is shown to be the non-linear ODE yy"=3x-(y')^2.

SOLUTIONS OF DIOPHANTINE EQUATIONS: Equations of the form a f(x)+b g(y)=c are referred to as Diophantine equations named after Diophantus of Alexandria who worked in Egypt during the third century AD. One is interested in finding only the integer values of x and y for which the equation is satisfied. After looking at the simplest case of linear Diophantine equations were f(x)=x and g(y)=y, we discuss non-linerar types such as the elliptic equations and the famous equation y^2-Kx^2=1 of Brahmagupta. Click on DIAPHANTUS to see the discusion.

SOLUTION OF THE NON-LINEAR ALGEBRAIC EQUATION Y

EVALUATION OF THE RIEMANN ZETA FUNCTION: We evaluate the Riemann Zeta function Zeta(s)=sum[1/n

RIGHT TRIANGLES AND THE BASE AND PYTHAGOREAN TRIPLETS: We examine the solutions to a

THE BABYLON PROBLEM OF IRREGULAR POLYGONS INSCRIBED WITHIN A CIRCLE: One of the oldest recorded geometrical problems in the history of mathematics is that of the properties of irregular N sided polygons confined to a circle and having all of its vertexes lying on the circle. This problem was first looked at by the Babylonians in about 2000BC and their cuneiform recorded tablets clearly show a knowledge of , among other things, the Pythagorean Theorem. It's really quite amazing how far they were able to progress using their rather clumsy sexagesimal number system with no understanding of zero. Go to BABYLON to see some variations on this N sided polygon problem.

DEFINING CIRCLES AND ELLIPSES BY THREE AND FOUR POINTS IN A PLANE: We show that three points in a plane, not lying along the same straight line, define a unique circle while four points are necessary to define an ellipse. An extention to 3D shows that four points not all lying in the same plane are required to define a sphere. Go to CURVES-DEFINED-BY-POINTS for the details.

CONSTRUCTING SYMMETRIC FIGURES USING DIFFERENT SUB-ELEMENTS: We show at N-FOLD-SUMMETRY how one can use rotations of simple sub-elements to generate m fold symmetric figures. Our computer calculations include the construction of a Greek upper case Delta, a Maltese Cross and an Eight Petal Flower. A computer generated Christmas greeting is also shown.

BESSEL FUNCTION IDENTITIES: The solution of numerous mathematical problems expressed in cylindrical coordinates involve Bessel Functions as part of their answer. At BESSEL-IDENTITIES we look at some of the more important identities involving these functions. It is shown, for example, that J

ITERATION TECHNIQUES FOR FINDING ZEROS OF FUNCTIONS AND POWERS OF NUMBERS: Most hand calculators and electronic computers have built into them programs which find zeros of functions f(x) and the value of numbers N taken to their pth power. This is accomplished by iteration programs using a Tayor Series as their starting point. We demonstrate at ZEROS-OF-FUNCTIONS how this works for several different cases.

2011

GEOMETRIC PROPERTIES OF PENTAGRAMS: We derive various sub-lengths of a standard pentagram. Go to PENTAGRAMS to see the details. The calculations involve use of the Golden Ratio. It is also shown how some of the triangular sub-elements may be used to produce aperiodic Penrose Tilings. Areas of the circumscibed and inscribed pentagons are also found.

PATHS TO CONVERGENCE FOR THE MANDELBROT SET: It is known that the Mandelbrot iteration A[n+1]=(A[n])^2+c, where c is a complex number leads to interesting color patterns when plotting different values of the function W=exp(-abs(A[inf])). What has not been considered in detail is the paths to convergence the A[n]s take when the iteration is bounded. We discuss at MANDELBROTthe type of paths to convergence in the A[n]=Re(A[n])+I Im(A[n]) plane as c=x+Iy is varied. It is shown that in most cases the paths to a convergence point A[inf] are spirals. Point plots are presented for some of the more interesting cases near the convergence-divergence boundaries.

SEVEN BRIDGES OF KOENIGSBERG AND RELATED ROBLEMS IN GRAPH THEORY: A famous problem in Graph Theory deals with whether or not it is possible to cross all the seven bridges across the river Pregel in Koenigsberg, East Prussia and end up at the starting point without crossing any bridge twice. Euler in 1735 was the first person to give a correct mathematical answer to this problem by showing that is not possible. In the process he invented Graph Theory. Go to KOENIGSBERG-BRIDGES to see a discussion of this problem and some other problems related to it.

INTEGER SPIRAL: Several years ago we noticed that the Archimedes Spiral r=(4/Pi)theta has integer values n at angle theta=(Pi/4)n and that r=n there. Thus the number 37 is located at r=37 at an angle of theta=37 Pi/4. All odd integers are shown to lie along the diagonal lines y=x and y=-x while even numbers are found along the x or y axis. Go to INTEGER-SPIRAL for details in which we show the location of the first 47 integers along the Integer Spiral and point out the location of prime and composite numbers. A discussion on addition and multiplication of integers using modular arithmetic based on mod 8 is also presented.

GENERALIZATION OF THE FIBONACCI SEQUENCE: Go to FIBONACCI-GENERALIZED to see our discussion on the Fibonacci Sequence and introduce a variation in which the g

PROPERTIES AND SUMS OF INTEGERS AND OF PRIMES: We look at the integer sequence 1,2 3,4 ,5,6,7,8,9,.. and sum the even , odd, and prime integers taken from this sequence. The sum of the first N primes is given approximately by N

`,`

n=1..infinity],
where
pMORE ON ARCTAN APPROXIMATIONS: We use the integrals int(P(2n,x)/(x^2+a^2),x=0..1)=N(n,a)+(M(n,a)/a)arctan(1/a) and int(P(2n+1,x)/(x^2+a^2),x=0..1)=U(n,a)+V(n,a)ln[1+(1/a^2)]

involving Legendre polynomials P(n,x) to obtain some excellent approximations for arctan(1/a) and ln(1+a^2)-ln(a^2) for larger a. Estimates for the accuracy of such approximations are given and specific values for arctan(1/5), arctan(1/38) and ln(2) are presented. Go to AGAIN-ARCTAN to see the discussion.

SYMMETRIC PATTERNS PRODUCED BY THE ITERATION Z[n+1]=Sum{C

KOCH CURVES: We show how a basic Koch curve is constructed and why it has infinite length but finite area between itself and the x axis. Some area Koch patterns are also generated using a Swiss Cross and a simple square as the starting point for the iteration leading to the final boundary. Go to KOCH-CURVESto see the details.

GENERATING LARGE PRIMES FOR PUBLIC KEYS: It is well known that public keys in RSA cryptography rely on the product of two large primes whose product cannot be factored with even the fastest existing supercomputers. Typically such secure public keys require prime nber components in the 1000 bit range and these are typically obtained by a random number search. Go to GENERATING-PUBLIC-KEYS where we discuss an alternative and faster prime number generation technique based on a linear combination of integers taken to a high powers and made prime by varying a parameter A. The approch is related to the Mersenne Prime formula but is more complicated and hence making it nearly impossible to guess which prime numbers are involved. Three large primes( among many others) found are N

FERMAT'S LITTLE THEOREM AND NUMBER PRIMALITY: A number is prime if it can be devided only by itself and by one. The standard approch to testing for the primality of a number involves variations on Fermat's observation that a number p is prime if the quotient [a^(p-1)-1)]/p has integer value. Typically a is taken as 2 but can have other values such as 3, 4, etc. At LTTLE-THEOREM-FERMAT we examine Fermat's Little Theorem and also look at a variation which states that p is prime if Sum[ k

FACTORING N=pq: It is well known that RSA cryptography depends on the inability for an adversary to factor large composite public key number N into the product of two primes p and q contained in the sender's private key. The present procedures for factoring large N involves essentially the brute force method of dividing N by all primes from 3 through the nearest prime to sqrt(N) or using a generalized numerical sieve or elliptic curve factorization procedure. All these approaches can take a very long time even with the fastest electronic computers. We look at FACTORING N to see the first and most primitive factoring procedure and show how one can accelerate somewhat this brute force approach by breaking the unknown primes into a form q=K-beta and p=K+alpha for a chosen K and then solving a resultant Diophantine Equation. The approach is a little faster than the brute force approach, but is still insufficient when talking about factoring very large Ns such as the example given at the end of the discussion.

FINDING LARGE PRIMES: For those cryptograpic methods where one employees a public key one is very much interested in using the product of two large primes P and Q. These primes should be about 200 to 2000 digits long each and not be too close in value to each other. We show here that a modification of the classic Mersenne formula may be used to quickly generate a large number of such primes . It is shown in detail how one can generate a 383 digit long prime in just a few minutes and then store or transmit it in the understandible compact form P=P[-14,(1,1,1),(7,13,17),(97,63,311)]. The generating formula used is P=A+sum[c

OBTAINING POLYNOMIAL QUOTIENT APPROXIMATIONS TO FUNCTIONS USING LEGENDRE POLYNOMIALS: We examine integrals of the type Int[P(2n,x)F(a,x),x=0..1] where P(2n,x) are high order LegendrePolynomkials and F(a,x) is a slowly varying function of 'a'. We find excellent approximations for the exponential function exp(a), the arctan function arctan(1/a), and the log function ln[1+(1/a)]. Go to POLY-QUOTIENTS to see the details of the procedure and a discussion of its limits. Among several other things, we show that

SOME EXACT VALUES FOR ARCTAN(1/a): We employ some known formulas for the arctan function to determine the exact values of arctan(1/a) for several different values of a. These calulations include a demonstrationt that Pi/16 can be represented precisely by arctan[(a-1)/(a+1)] when a={1+sqrt[1+(1+sqrt(2)^2]}/[1+sqrt(2)]. Also the expansion of arctan(1/a) as a finite series of other arctan functions with larger values of 'a' are derived. Finally we look at the arctangent of the complex variable z=x+iy. A contourplot of arctan(z) shows a pattern reminiscent of a source-sinhk flow in fluid dynamics. Go to EXACT-ARCTAN for the details.

MORE ON FUNCTION APPROXIMATIONS USING LEGENDRE POLYNOMIALS: We have shown earlier that arctan(1/a) can be well approximated by a polynomial quotient in a with accuracy increasing as n in the even Legendre polynomials P(2n,x) used in the method is increased. At FUNC-APPROX we extend this discussion to obtain approximations for more functions including exp(2a), sinh(a), cosh(a), tanh(a), erf(a) and Si(a).

CONSTRUCTING CLOSED CHAINS OF CIRCLES: A computer program is employeed to construct closed chain configurations of circles and moment of inertia concepts are then applied to measure their compactness. Go to CHAINS-CIRCLES for the details. It is shown that such rings of constant radius circles is achieved by anchoring their centers to the vertexes of regular n sided polygons and that their radii R are half the polygon side-length of Ssin(Pi/n), where S is the distance from a circle center to the origin. We also show that the largest circle which may be placed in the center of the circular ring chain region is S[1-sin(Pi/n)].

EVALUATING TAN(a) TO ANY ORDER OF ACCURACY: In a recent note we showed how one can approximate the values of the trignometric functions using an integral approch involving the odd Legendre polynomials. We extend this discussion at HIGHLY-ACCURATE-TAN by showing how one can achieve very high order accurate approximations for the tangent function tan(a) over the entire range -inf<a<+inf.

Also the location of its singularities at (2k+1)Pi can be accurately determined. A fifth order approximation involving the quotient of an eleventh order polynomial divided by a tenth order polynomial in 'a' produce forty digit or greater accurate values for tan(a).

SOLUTION OF THE BRAHMAGUPTA-PELL EQUATION: We examine the equation y

PROPERTIES OF THE SINE INTEGRAL: We examine the area under the curve sin(x)/x extending from 0 to x This involves the sine integral Si(x). After giving its usual series expansion we look at approximations for both large and small x. Also an integral representation for the Si(x) function is derived and then evaluated using an approximation for the tan(t) function. Go to PROPERTIES-SINE-INTEGRALfor the details.

SUMMING THE RECIPROCAL OF THE FIRST N INTEGERS: Unlike summing the sum of the first N integers, the task of finding S(N)=Sum[1/n,n=1..N] is a bit more complicated and involves the Digamma function Psi(N). We show at RECPROCAL-INTEGERS how this sum is obtained and prove that the sum equals S(N)=Psi(N+1)-Psi(1). In the development we find the interesting identity that

Int[(1/x-1/t)t^x ln(t) exp(-t), t=0..infinity]=(x-1)!/x. Using Laplace transforms we find the value of Psi(1) to equal the negative of the Euler constant gamma=0.577215.. . All other values Psi(x) can be obtained from the recurrence relation Psi(x+1)-Psi(x)=1/x.

EVALUATION OF FUNCTIONS BY ITERATIONS BASED ON THEIR SERIES: I have often wondered exactly what procedures my hand calculator and PC math programs use to find values of a function at any given point in a split second. Clearly the procedure these built-in programs are using is some form of iteration based on the series representation of the functions. We show at FUNCTIONS-BY-ITERATION how this is probably done . Specifically we look at iterations for values of exp(x), arctan(x), and Si(x). Furthermore we show that a very rapidly converging iteration for sqrt(2) is given by S[n+1]=A/B-1/(B(A+BS[n])). The constants have the integer values A=886731088897 and B=627013566048.

SUMMARY OF OUR WORK AN ARCTAN AND TAN APPROXIMATIONS: During the summer and fall of 2011 I collaborated with Sidey P. Timmins of Sabre Systems at the US Census Bureau to summarize, clarify, an extend some of the work I have done recently on the use of integrals involving Legendre Polynomials to accurately approximate the arctan and other trignometric functions by Pade like quotients. He has been particularly interested in incorporating the procedure into standard computer routines requiring orders of accuracy in excess of 40 decimal places. Attached you will find two preliminary papers we have written on the topic. The first (APPROX-ARCTANH) deals with arctan approximations and the second (APPROX-TAN) with approximations to the tan function and related trignometric functions.

FINDING THE SQUARE ROOT OF NUMBERS BY AN ACCELERATED ITERATION TECHNIQUE: We look at the problem of finding the square root of numbers and discuss how these were obtained in pre-computer days. Also , we introduce an accelerated iteration technique based on continuous fractions and the fact that [sqrt(N)-sqrt(N

THE HOUGH TRANSFORM: In the early 1960s Paul Hough introduced a digital analysis technique which is capable of finding all points in a 2D array which fall along the same straight line. Now known as the Hough Transform, it is one of the main methods used in digital image processing. We discuss at HOUGH-TRANSFORM how this transform works and give an example of finding the line along which four pre-specified points lie. Also we consider an even more elementary geometrical method for finding straight lines connecting more than two points within an array of N points.

SUM AND PRODUCTS OF PRIMES: We examine the sums and products of the primes 2, 3, 5, 7, 11, 13,. ..It is shown that sums of the form F(s,x)=Sum[1/p

POWERS OF REAL AND COMPLEX NUMBERS: We evaluate numbers of the form N

2012

EVALUATING EXP(X) AND RELATED FUNCTIONS TO 50+ DIGIT ACCURACIES: In several earlier discussions we have shown how to generate good approximations to the values of different trignometric functions using a new technique based on integrals involving the standard Legendre polynomials. We extend at 50-DIGIT-ACURACY these investigations by showing how highly accurate approximations to exp(x) and variations thereof can be obtained by using an integral involving the product of P(2n,x) and cosh(x/2a) integrated over the range [0,1].

EPICYCLOIDS, HYPOCYCLOIDS, EPITROCHOIDS, AND HYPOTROCHOIDS: An interesting set of curves can be generated by roling a circle about the outside or inside of a stationary circle. Go to EPICYCLOIDS to see how one arrives at parametric equations for the four distinct curves resulting when tracing out the path of a point P which lies along a radial line of the rolling circle at fixed distance from the circle center. Graphs of the better known configurations are presented including the Astroid.

INTEGER SPIRAL: Several years ago we first noticed that any positive integer can be located along the Archimedes Spiral r=(4/Pi)theta at the intersections with the x or y axis or with the diagonals x=y or x=-y. We call this spiral the Integer Spiral and have shown earlier that the well known Ulam Spiral easily morphs into the present spiral form thus negating any thoughts that the Ulam Spiral contains hiden information involving the prime numbers. Modular arithmetic is used to locate any integer in the x-y plane. Also, we show that all Mersenne Numbers lie along a diagonal in the 4th quadrant of the x-y plane while all Fermat numbers lie along the diagonal in the first quadrant. Go to INTEGER-SPIRAL for the details of the discussions.

GOLDBACH CONJECTURE: In 1742 Christian Goldbach proposed that twice the product of any positive integer N is equal to the sum of two primes Pn and Pm. Known as the Goldbach Conjecture it remains one of the most important and unproven theories in number theory. We discuss at GOLDBACH how to find a pair of primes whose sum equals 2N, and show how one can use Goldbach's conjecture to quickly break a public key R=PnPm by varying N near sqrt(R).

PROPERTIES OF THE NUMBER TRIANGLE: It one rotates a standard multiplication table by 45deg there will result a number triangle of the form-

1

2 2

3 4 3

4 6 6 4

5 8 9 8 5

6 10 12 12 10 6

Go to TRIANGLE to see the details of the properties of such a triangle and show that K(n,m)=m(n+1-m) yields the value of all its elements with n being the row number and m the column number. The value of all elements along a given row fall on a unique parabola. The sum S(n) of these elements in a given row n is given by n(n+1)(n+2)/6. Each odd number within the Number Triangle is found to be surrounded by six even numbers which are easily determined once the central odd number is specified.

CONSTRUCTION OF SPIRALS USING STRAIGHT LINE SEGMENTS: We examine the type of inward and outward winding spirals which can be constructed by connecting straight line segments to each other. Go to SPIRAL-CONSTRUCTION for the details. Among other spirals we generate the Number Spiral and discuss some of its properties. Other spirals are also generated including the Ulam Spiral and a spiral which approaches the Archimedes spiral as the angle increment in a polar plot becomes small. In addition, we look at a few cases where the radial distance from the origin to the line segment is periodic. This produces some interesting non-spiral symmetric curves including those associated with fractals.

GENETIC ALGORITHMS FOR GENERATING 2D CURVES: Just as a DNA molecule stores genetic information as a sequential arrangement of four bases, we show you at GENETICS that the same thing can be done with generating 2D mathematical curves using unit length straight line segments as the building blocks subjected to just three bases labelled Left(L), Right(R), and NoTurn(O). For example any square can be defined by a four time repetition of the single letter code- R-. The code O-R-R-L-O-R-R-O-L-R-R-O produces the closed curve representing the capital letter T. Note that this last code has the form of a palindrome predicting that the letter T has a symmetry line along its vertical axis. In addition to right angle bends between straight line segments we also examine the posibility of any other angle bends as encountered with regular polygons such as equilateral triangles or hexagons.

EVALUATION OF A GENERALIZED GEOMETRIC SERIES: We examine the generalized geometric series f(z)+f(z)

CONSTRUCTING N POINTED STARS: The construction of N pointed stars using the listplot option in our MAPLE program is demonstrated. Symmetric stars with five through ten points are created and formulas given for extending these calculations to any N pointed star. In addition we look at several other configurations including a twenty pointed star. Go to POINTED-STARS for the details.

USE OF THE GOLDBACH CONJECTURE FOR FACTORING LARGE SEMI-PRIMES: It is known that it is very difficult to factor large digit composite numbers N which equal the product of just two primes p and q. We take another look at this problem at GOLDBACH-AND-SEMIPRIMES and introduce the Goldbach Conjecture to aid in our efforts. We find that p=n-sqrt(n^2-N) and q=n+sqrt(n^2-N) , where n is the mean value of the unknowns p and q. Treating n as a running variable we show how a large digit semi-prime N can often be factored into its two components with little effort. We also factor the Fermat number N=2^32+1 by making use of the fact that N/n^2 is a very small number in this case.

GENERATING PERFECT SQUARES FROM THE PRODUCTS OF NON-SQUARES: A very important manipulation used in the factoring of numbers is to construct perfect squares from a collection of non-square integers. One way to do this is to add together the elements of various exponent vectors based on powers encountered in prime number factorization. As an example consider the product of the non-square numbers 40 and 90 . It produces the perfect square 3600=(60)

THE GREATEST COMMON DIVISOR: The operation gcd(N,M) refers to the greatest common divisor between two integers N and M. It represents essentially the largest divisor which factors into both N and M exactly. For example, gcd(72,45)=9 since 9 is the largest number which divides both 72=8x9 and 45-5x9. We discuss at GCD the several techniques available for quickly finding gcds and then apply the operation to solve the linear Diophantine Equation and to factor a semi-prime number.

AN ALTERNATIVE METHOD FOR FACTORING LARGE SEMI-PRIMES: We show how a semi-prime N=pq can be factored into its two components p and q by the evaluation of the new formula gcd((n

EVALUATION OF THE COMPLETE ELLIPTIC INTEGRALS BY THE AGM METHOD: We examine the values of the complete Elliptic Integrals K(m) and E(m) using Gauss's Arithmetic-Geometric Mean Method(AGM). After discussing the variable transformations which make the AGM approach possible, we show details of the required iteration approach. Go to ELLIPTIC-INTEGRALS-BY=AGMfor the details. Among other points we obtain the value of the perimeter of a lemniscate accurate to 87 decimal places.

EVALUATION OF CERTAIN DEFINITE INTEGRALS RELATED TO THOSE SOLVABLE BY AGM METHODS: We examine certain intergrals of the form int[1/((x^2+1)(x^2+4)(x^2+9)...)), x=0..infinity. Closed form solutions are found . All of these solutions are shown to be directly proportional to the first power of Pi but unlike the integral int(1/sqrt((a^2+x^2)(b^2+x^2)),x=0..infinity) do not lend themselves for solutions by the AGM method. Go to INTEGRALS-VIA-AGM for the details of the discussions including the demonstration that the integral int(1/(x^2+1)(x^2+4)(.....)(x^2+64),x=0..infinity)

=Pi/6096384000.

ANOTHER ARCTAN FORMULA FOR PI: Although numerous arctan formulas which generate the value of Pi exist, we look here at a new single term arctan formula based on the equality Pi=(4 or 6)2

MERSENNE NUMBERS AND THE INTEGER SPIRAL: We examine the Mersenne Numbers M[n]=2

GRAPHICAL INTERPRETATION OF THE SIEVE OF ERATOSTHENES: One of the simplest ways to generate low number primes is the Sieve Method of Eratosthenes which starts with a list of the ascending odd numbers and then eliminates those numbers divisible by the first n prime numbers. This procedure is equivalent to looking for the zeros of the sine function sin[(x Pi /ithprime(n)] for a set of the lowest prime numbers n=3,5,7,11,13,..,229. Go to SIEVE for details of the discussion in which we use the sine curves of the first fifty primes to find prime numbers up to eleven digits in length by a simple visible inspection of the sine graphs. We show, for example, that 71423746921 and 89645676767 are prime numbers.

GENERATING INFINITE PRODUCTS OF FUNCTIONS BY EULER'S METHOD: We show how one can represent evenly symmetric functions having an infinite number of zeros as infinite products following essentially the approach first used by Leonard Euler. In particular we concentrate on the functions sin(x)/x, cos(x), tan(x)/x, and the Bessel Function J

MORE ON THE SOLUTION OF THE BRAHMAGUPTA EQUATION: The non-linear Diophantine equation y

PROPERTIES OF THE INTEGER SPIRAL: We look further into the properties of the Integer Spiral with which positve integers can be represented as points of intersection of an Archimedes Spiral and straight lines emanating from the origin. Using mod operations we locate numbers such as Mersenne Numbers and Fermat Numbers in the Argand Plane and show how they can be factored. Go to INTEGER-SPIRAL-PROPERTIES for the details.

CONSTRUCTION OF FIBONACCI SPIRALS: In the classic construction of a Fibonacci Spiral one places different size squares tangent to each other and then draws circular arcs from one corner of each of the squares. The result is a Fibonacci Spiral which is continuous but lacks continuity in its 2nd and higher derivatives. We show at FIBONACCI-SPIRAL how to correct this shortcoming by using an approximation to the Binet Formula for large n to find the continuous Fibonacci Spiral r=0.723606 exp(0.6126979 theta ). The radial distance from the origin to this curve for a given n very closely matches the values of the Fibonacci Numbers F[n].

NUMBER FRACTION FOR INTEGERS: We examine the integers and introduce a new parameter termed the number fraction f(N) to distinguish composite from prime numbers. This number fraction is defined as-

where sigma(N) is the sigma function of number theory representing the sum of all divisors of N. Prime numbers correspond to f(N)=0. It is also shown that numbers Q=6(random number+n)$\pm $1 , where n represents certain integer values, will be primes. Also we find an important sub-class of primes given by Q=6 x 2

FURTHER DISCUSSIONS ON GENETIC CODES FOR 2D CURVES: In an earlier note we showed how certain 2D curves can be generated by simple building blocks not unlike what happens in genetic DNA reproduction. We extend this discussion at MORE-GENETIC-GENERATED-CURVES to base elements of variable length and variable angles between the elements. Very elaborate curves are shown to follow from simple base elements. Among the figures generated is a square spiral, a five pointed star, and a periodic pulse function. Also it is shown how such straight line elements can produce continuous curves under appropriate limiting conditions.

GENERATING Q PRIMES: In a recent note we discussed how to distinguish between prime and composite numbers using the concept of the integer fraction. In that study we observed that integers with the largest number of factors are those divisible by 6 and 12 . Furthermore we discovered that all primes except 2 and 3 are given by the prime versions of Q=6N$\pm $1. We examine these primes in more detail a t HGENERATING-Q-PRIMES . The first hundred Q primes are listed and then we demonstrate how such Q primes of hundred of digit length are generated with very little effort. Also it is shown how to break large public keys N=pq when both p and q are Q primes. The latter may find application in cryptography since the Q primes constitute essentially all primes above 3.

MORE ON FACTORING LARGE PRIMES: In several earlier articles we have discussed methods for factoring large semi-prime numbers. Go tot FACTORING-LARGE-PRIMES for discussions and seeh , among other things, why the factoring of N=pq involves essentially a solution of the Brahmagupta-Pell equation . By denoting M as the average value of p and q one needs to search for that value of M which makes d=sqrt(M^2-N) equal to an integer. Once this has been done, one has the solution p=M+d and q=M-d. It is shown how the eleven digit number N=79998475553 is factored into its prime number components p=298621 and q=267893.

USING THE INTEGER SPIRAL TO FACTOR SEMI-PRIMES: We want in this note to show how one can locate semi-primes N and their factors p and q along diagonals in the integer spiral and use the result to find p and q knowing only the quantities ab=N mod(8) and K=(N-ab)/8 in the expression N=(8n-a)(8n-b). The problem boils down to finding the value of V=am+bn for which sqrt[4V^2-2ab(V+K)] is an integer. Go to INTEGER-SPIRALS-AND-PRIMES for the details.

RELATIONSHIP BETWEEN PERFECT NUMBERS AND THE NUMBER FRACTION: We discuss Perfect Numbers such as R=6, 28, and 496 and show that they are related to our recently defined Number Fraction f

PRIME NUMBER TEST: We show that all prime numbers N above N=3 must be of the form N mod(6)=1 or 5. Using this information together with an evaluation of the number fraction f

THE BINOMIAL EXPANSION AND ITS MULTIPLE VARIATIONS: We look at the classic Binomial Expansion (a+b)

APPROXIMATIONS FOR FUNCTIONS USING ITERATIVE APPROACHES: We show at ITERATIVE-APPROACH how one can use iteration methods based on the Taylor Series to approximate the values of functions such as x

DETERMINING VALUES OF ARCTAN(1/N) BY ITERATION- We extend our discussions on iteration methods by looking in more detail at arctan(x) functions, where N=1/x is a large number. It is shown that less than ten iterations will often produce values for arctan(1/N) which are accurate to well over 80 digits. The iterative algorithm for doing so reads z[n+1]=z[n]+(1/N)cos(z[n])

FERMAT-PYTHAGOREAN-TRIPLETS: We look at the integer solutions to a modified Fermat-Pythagorean Formula a

CONSTRUCTION OF 2D OPEN AND CLOSED CURVES USING GENETIC CODES: We examine the properties of curves consisting of unit length segments connected to each other at fixed angles of +pi/2, 0, or -pi/2 and represented by genetic codes such as [+ - 0 - +]. It is shown that when the number of + signs in the code equals the number of - sings, the curve generated by n application will generally produce an open curve. When the number of + and - signs differ from each other, then the code will likely produce a closed curve upon multiple applications, although exceptions do occur. We show that four concatenations of the genetic code [ - + + ] produce a Swiss Cross. The code [ + 0 + ] will produce a rectangle of two to one aspect ratio. Go to CURVES-BY-GENETIC-CODESfor the details.

2013

SOLUTION OF CERTAIN NON-LINEAR ALGEBRAIC EQUATIONS: We examine the non-linear algebraic equations governing the length of the hypotenuse (longest side) of an isosceles triangle formed by connecting the nth and n+2 corner of a regular n sided polygon. It is shown that these highly non-linear equations will always have 2cos(Pi/n) as one of its n solutions. Go to NON-LINEAR-EQS for the details of the discussion. It is also found that sin(mx)/sin(x)=sum[(-1)^2(m-k-1)!(2cos(x))^(m-2k-1)/(k!(m-2k-1)!), k=0..(m-1)/2] for m=2, 3, 4, etc. . A closed form expansion of cos(mx) in terms of powers of cosine is also obtained.

** SEMI-PRIMES AND THE NUMBER
FRACTION: The number fraction for a semi-prime N=pq
is given by f=(p+q)/N and can be easily generated by a PC for
Ns up to about 50 digit length. We look at three specific
semi-primes, find their f values when possible, and then give
explicit values for the prime number factors p and q. We
show, for example, that-N=**

**GENERATING PERFECT NUMBERS:
We examine perfect numbers N defined as
b=Sum(all divisors)-2N=0. The lowest few of
these are N=6, 28, 496, 8128, 33550336.They
relate to Mersenne Numbers M =2 ^{p}-1
via the identity N=M x 2^{(p-1)} where p
is a prime for which M is a prime. We also show
how N and M lie along separate radial lines when
plotted as part of an integer spiral. Go to PERFECT-NUMBERS
for the details. **

SEMI-PRIMES AND THE
SOLUTION OF Y=SQRT(X^{2}-N): It
is shown how the semi-prime N=pq is equivalent
to the Diophantine equation y^{2}=x^{2}-N.
A graph of this equation is presented and it is
then shown how one goes about finding the
desired integher solution pair [x_{1,} y_{1}]
=[(p+q)/2, (p-q)/2]. Specifically we show
how N=12063943 factors into p=5171 and
q=2333. Go to SOLVING Y^2-X^2-N for
the details of the discussion.

AN ALTERNATIVE METHOD FOR
FACTORING SEMI-PRIMES: We examine a
new approach for factoring semi-primes N=pq
based on the tangent of an angle in a
sqrt(N)--(p-q)/2--(p+q)/2 right triangle. It is
shown that p is a factor of N only if the
quantity F(N,p)=|(p^{2}-N)/(2p)| has
integer value. Thus N(119,7)=5 but
N(119,13)=25/13, meaning p=7 is a prime factor
of 119 but p=13 cannot be. Plotting F(N,p)
versus p produces a two branch curve which
merges to zero at p=sqrt(N). Semi-primes of up
to ten digit length are easily factored by this
approach. Go to
FACTORING-N=pq for
the details.

MORE ON SOLVING
F(N,x)=|(x^2-N)/(2x)| :
We conclude our discussions on factoring
semi-primes N=pq held over the last several
years by looking at solving the Diophantine
equation F(N,x)=(N-x^2)/(2x) for integer x
and F for fixed N. Such solutions produce
the integer factors q and p for any
semi-prime. We study in some detail the
F(N,x) graph associated with the semi-prime
N=2047 whose factors are q=23 and p=89. Also
we indicate how stored values for the first
50,000 primes can be used to quickly solve
F(N,x) as long as N is less than about 12
digits in length. Go to
SOLVING-A-QUOTIENT for the details.

N! AND THE GAMMA FUNCTION: We look more at the factorial and gamma functions and discuss many of their properties

and present graphs whereΓis a function of the complex variable z=x+iy. Go to N!-AND-GAMMA for the details.The binomial coefficient , the

Psi function, sums of reciprocal powers of n! , and the function P=Γ(x+y)Γ(x-y) are examined in some detail. A general formula for

the value of all half-integer Gamma Functions is derived based on the Legendre Duplication Formula.

DETERMINING WHETHER A NUMBER IS PRIME: It is shown that any prime number above N=3 can be represented as

Q=6n+1 or 6n-1where n=1,2,3... Prime numbers of 5 or greater have the propertythat N mod(6)=1 or 5, howeverfor n=1,2,3..

this is not sufficient unless the ratio test R=N/(6n+1) or N/(6n-1)fails to yield a integer value. Any odd number for which N mod(6)=3 will always be a composite.Thus N=456723108745338992349075615353107312371is a composite. Go to

DETERMINING-IF-N-IS-PRIME for the details of the procedure for factoring large composite numbers into products

of primes.

PRIMES, COMPOSITE, AND PERFECT NUMBERS:We use the Number Fraction f(N)=[sigma(N)-N-1]/N to distinguish

between prime and composite numbers and show that f(N)=1-(1/N) produces the perfect numbers N=6, 28, 496, etc. Here

sigma(N) is the standard sigma function defined as the sum of all the divisors of the number N including 1 and N. In particular

we look at super-composites which are numbers for which f(N)>1. The largest value of f(N) found in the range 1<N<700000

is f(665280) =3.39826. Go to PERFECT-NUMBERS for the details including the facts that f(p^{n})=[1-(1/p^{n-1})]/(p-1), where p is any prime number,

and that the two numbers in the immediate vicinity of a super-composite tend often to be prime numbers.

**MORE ON
FACTORING SEMI-PRIMES
USING NUMBER FRACTIONS: We
examine the number
fraction of semi-primes
in the neighborhood of
large supercomposites.It
is found that such
semi-primes x=pq can be
factored into their
components by p=a+sqrt(a ^{2}-x),
where a=xf/2, with
f=[sigma(x)-x-1]/x being
the number fraction
andsigma(x) the sigma
function. Three specific
semi-primes , namely,
4047, 2^32+1, and
1540(6^49)-1 are treated
in detail. In particular
we show how the
factoring-**

**PROPERTIES OF N SIDED REGULAR POLYGONS: We
examine the
general
problem of
regular convex
polygons and
show how many
can be
constructed
using only
straight edge
and compass.
All vertex
points of such
polygons lie
on a circle of
constant
radius with
neighboring
points
separated from
each other by
a vertex angle
of 2Pi/n, with
n being the
numbers of
edges of the
polygon. The
17 sided
polygon of
Gauss is
discussed and
we give a
detailed
description
for the
construction
of a regular
pentagon. It
is also shown
how any
regular
polygon can be
readily
magnified,
translated and
rotated in the
x-y plane. Go
to
N-SIDED-POLY
for the
details.
**

**FINDING THE CIRCLE WHICH INSCRIBES ANY TRIANGLE:
We find the
center [a,b]
and radius r
of a circle
which passes
through all
three vertices
of a given
triangle. The
procedure
involves some
simple matrix
algebra
involving the
inversion of a
2 x 2
coefficient
matrix M and
works as long
as det M does
not vanish. Go
to CIRCUMSCRIBING-CIRCLE for
the details.
It is also
shown how a
magnification
and rotation
of a given
triangle can
be quickly
accomplished
by using
information on
the
circumscribed
circle.**

ROTATION
OF ANY 2D
CURVE: We
examine the
problem of
rotating a
general two
dimensional
curve
F(x,y)=ax^{2}+2xy+cy^{2}
+ex
+fy+h=0
by an angle
theta using
matrix
methods.
Defining
coefficient
matrixes C and
G and a
rotation
matrix M, we
show that the
2D curve
converts to a
new form X'^{T}[MCM]X'+GMX'=-h.
HereX' is the
column matrix
containing the
new variables
x' and y' and
the
superscript T
is the
transpose.
Several
different
examples are
looked at
including a
rotated
straight line,
a rotated
ellipse and
hyperbola,
plus several
rotations of a
standard
parabola. Go
to CURVE-ROTATION for the details.

EULER
ANGLES AND 3D
ROTATIONS BY
MATRIX
METHODS: We examine the procedure for quickly
performing
several
rotations on
objects using
3D Rotation
Matrixes.
After defining
the Euler
Angles, we
show how one
can take a
normal cube
and rotate it
two times to
produce a
orientation
where the cube
diagonal
coincides with
the vertical z
axis. Also we
show how to
transform any
point P(a,b,c]
to a new point
P[e,f,g] by
three
independent
rotations
about the x,
y, and z axis.
In addition,
we
demonstrate
how the
slanted plane
x+y+z=1 can be
rotated into a
new plane
z=1/sqrt(3) by
just two
rotations. Go
to EULER-ANGLES for the details.

MORE
ON Q PRIMES: We
discuss in
greater detail
when numbers
of the form
N=6n±1 are Q
Primes. A
simple test
consisting of
dividing N by
6k±1 and then
looking at the
values of this
quotient over
the integer
range
1<k<sqrt(N)/6
is sufficient
to determine
whether N is a
prime or
composite. If
no integer
values are
found then the
number N is a
Q Prime. If
one or more
integer values
occur than N
is a
composite. All
odd
numbers
having N mod
(6) =3 must
always be a
composite. We
also present a
new hexagonal
integer curve
in which all Q
primes lie
along just two
radial lines
at angles **±**Pi/3.
The curves
looks like
this-

**Go
to Q-PRIMES for further details.**

FINDING
LARGE
SEMI-PRIMES:
In public key
cryptography
one deals with
large
semi-primes in
excess of 100
digit length.
Such numbers
are difficult
to factor and
hence
guarantee
security of
encrypted
communications.
This will no
longer be the
case in the
future when
high speed
supercomputers
will be able
to factor such
numbers
in
reasonable
amounts of
time. It is
our purpose
here to
identify when
a number N is
a semi-prime.
It is shown
that such
numbers and
their
components p
and q all lie
along just two
radial lines
intersecting a
hexagonal
integer
spiral. One
can identify
large
semi-primes by
looking at the
value of its
number
fraction
f={sigma(N)-(N+1)}/N.
If the value
lies near
f=2/sqrt(N) ,
but is greater
than zero,
then it is
likely to be a
semi-prime. We
show that
N=460969682477
is such a
semi-prime. Go
to LARGE-SEMI-PRIMES for the details.

SOLUTION
OF DIFFERENCE
EQUATIONS IN
THE COMPLEX
PLANE:
We examine the
difference
equation
Z[n+1]=F{Z[n]}+a+Ib,
where F
represents a
function of
Z[n].
Trajectories
of the
solution paths
Z[n] in the
complex Z
plane are
determined for
several
different
cases. The
paths starting
with Z[0]=Z_{0}
often lead to
interesting
spiral
patterns
ending at a
point
Z[infinity]
which can be
determined
analytically
by noting that
Z[n+1]=Z[n] as
n becomes
infinite for
converging
solutions. Go
to
DIFFERENCE-COMPLEX for
the details of
the
discussion.

GENERATING
ITERATION
FORMULAS:
We show how
one can
generate
iteration
formulas for
various
constants
including
exp(1), Pi,
ln(2), the
Golden Ratio,
Euler-Mascheroni
Constant, and
the mth root
of the number
N. It is
observed that
when the first
term in the
iteration lies
close to the
desired final
value then the
number of
iterations
required for a
good
approximation
will usually
be quite
small. Go to GENERATING-ITERATION-FORMULAS for the
details. When
a constant is
defined by an
infinite
series , the
n+1 iteration
relates to the
nth iteration
as simply the
n+1 term in
the infinite
series. Thus
the constant
ln(3/2) can be
evaluated by
the iteration

S[n+1]=S[n]+{2/(2n+3)}5^{-(2n+3) }subject to S[0]=2/5.

COTES
THEOREM:
In the early
17 hundreds
Roger Cotes of
Cambridge
University
came up with a
theorem which
states that
the product of
the distances
between N
equally spaced
points on a
unit circle
and an
interior
point
Q equals
1-x^{N}.