Web Page Started November 20, 1995
                                                         Latest Update April 19, 2018

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 a4,2=15 and  a4,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.


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 =2^n*Gamma(n+1/2)/sqrt(Pi) is compared to an evaluation of =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[cnfn(t), n=1..N], where fn(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 fm(t) and integrates over the range 0<t<T, this will lead to a set of non-linear algebraic equations for the cns which when solved give an analytic approximation for x(t). We demonstrate this approach for x"=x2 with x(0)=x(1)=0. Choosing just a single trial function f1(t)=c1sin(pt) , one finds c1= -(3/8)p3. 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.


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 yt+Ryx +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)zxx+2b(x,y)zxy+c(x,y)zyy=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 zxx+xzyy=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 [M2-1]jxx -jyy=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[M2-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 an=(2/L)*Int[U(x,0)*sin(n*Pi*x/L),{x,0,L}] and bn =(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 Utt=c^2[U xx+Uyy] 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 Ut(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
G(x,y)=sin(n*pi*x/a)*sin(m*pi*y/b) which satisfy all of the boundary conditions. These are multiplied by a spatially dependent term T(t)=Anm*sin(w *t)+B nm*sin(w*t), with the A's and B's determined from the initial conditions . Here w=p *c *sqrt[(n/a)^2+(m/b)^2] is the eigenfrequency. We show here the eigenmode corresponding to n=4 and m=5.

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 Utt=c^2*[Urr+(2/r)*Ur ]. 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 Ut(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 Vtt=c2Vxx+F(x,t) subject to say V(0,t)=V(a,t)=0  and V(x,0)=Vt(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
C"+(np/a)2C=D(t). Solving this by standard methods then yields C(t) which is substituted into the above sum for V(x,t). By clicking on the above title you can see a computer evaluation of the infinite sum for V(0.5,t) for the case of a string stretching between x=0 and x=1 when c=1 and F(x,t)=sin(wt) and there is no initial displacement or velocity on the string. Note the large response as one approaches the point where the driving frequency matches one of the natural frequencies.

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 Tt=Txx 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 Tbar=[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=aTxx 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 Tt= a [Trr +(2/r)Tr]  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 2t), this leads to R"+l 2R=0. From this follows the closed form solution T(r,t)=Sum[(2T oa/npr)(-1)n+1 sin(npr/a) exp( a(np/a)2t), 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 8cm/ p ) 2/0.205=2.01x10 17sec=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 Tt=aTxx+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)2C(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)2t)*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 +fyy=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 +fyy=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 fxx +fyy=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 F1(z)=Uoz and a source flow solution F2(z)=[Q/(2p)]ln(z). in this case the combined streamfunction is y(x,y)=Uoy+[G/(2p)] arctan(y/x). A contour graph for these streamlines when G/(2pUo)=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)=Uo(z+a2/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 zo=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[Anm(r/a)n*Ynm(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 Pn(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[Pn(x)Pm(x), x=-1..1]=2/(2n+1), are A0,0= 2/3 and
A 2,0= -2/3. That is, the temperature at any point inside the sphere equals T(r,q)=(2/3)-(1/3)r2( 3(cosq)2 -1). Click on the above title to see how isotherms corresponding to this solution look. As expected from the theorem of the mean, the temperature at the sphere center is 2/3.

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 ( x0,y 0) in y>0, and G is constructed from two line charges of opposite sign placed at P(x0 ,y 0)and P'(x0,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 x0=0, y 0 =1 is G=-[1/4*Pi)]*ln{[(x-x 0)^2+(y-y 0)^2]/[(x-x 0)^2+(y+y0 )^2]}. The solution of the Laplace equation thus becomes the Poisson integral f (x 0,y0)=[(1/Pi)]*Int[ f(x) dx/((x-x0 )^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 yxx+k2 y=0 with solutions yR=A Rexp(ik1x) and yI =A Iexp(-ik1x) for the reflected and incoming wave at the left of the interface. Also the solution is y T =ATexp(-ik2x) 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, yx ,  and the pressure , which by Bernoulli reads -i rwy, are continuous for the small amplitude waves considered.  Eliminating AT from the two resultant equations, yields the reflection coeficient R=(AT /AI)^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 c2)/( r 1c 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(xo,yo,z o)=G*Int[ r(x,y,z)dxdydz/r], where r=sqrt[(x-xo)^2+(y-y o )^2+(z-zo)^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/zo 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
f = [(3*M*G)/(2*a)]*[1-(1/3)*(z o/a)^2]. A plot of f over the entire range of zo is shown in the accompanying graph. Note the continuity of f at r=a and the zero slope of f at r=0 indicating that there is no gravitational force( F is proportional to grad f) at the sphere center.

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[Ro ^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/Ro )*[(r/R o)-2*cos( q )]. This will be a rapidly converging series whenever e<<1 and thus when the distance Ro 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/Ro , 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/Ro , 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 x1 and that from the cg to the moon center be x2, we choose a length scale where x 1 +x2=1 for the earth-moon distance. Now the total gravitational potential felt at a point P at distance ro from the cg is GM/rM+Gm/r m , where r M is the M-P distance and r 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 w2x 1 or the moon's centrifugal force  mw2x 2 , that  w2=GM/(x2(x 1 +x 2)2)
= Gm/(x1(x 1 +x 2)2). From this last result one also has that M*x1=m*x 2. Now the centrifugal potential at point P (x,y) due to a unit mass there is given by w2 ro 2/2 . Combining this last term with the two gravitational terms, one finds that the total potential felt by the test mass at P(x,y) is:

j(x,y)=GM{1/sqrt[x^2+y^2+x 1^2+2*x1*x]+(m/M)/sqrt[x^2+y^2+x2^2-2*x2 *y]+(x^2+y^2)/(2*x2*(x1+x2 )^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 yxxxx+2*yxxyy+yyyyy=0, with y being the streamfunction. Any solution of this equation represents a biharmonic function but it is only a special subclass of these functions which also will simultaneously satisfy an appropriate set of boundary conditions. There is a general theorem(see Sneddon's book on PDEs) which states that a biharmonic function can always be generated from y=xf1(x,y)+yf2(x,y)+f3(x,y)+f4(x,y), where fns are any harmonic functions satisfying the Laplace equation. Recalling how earlier we had a complex velocity potential F(z) for invicid 2D flows, it seems that one should be able to construct biharmonic functions by simply taking the real or imaginary part of F(z)=zbar*F(z)+z*G(z) +H(z)+J(z), where zbar is the complex conjugate of z=x+I*y and F(z), G(z), H(z) and J(z) are  analytic functions of a complex variable. To demonstrate this point consider the function constructed form F(z)=1/z, G(z)=1/z^4, H(z)=J(z)=0. That is F(z)=zbar/z+1/z^3. Taking the real part, we find the biharmonic function  y=[x(x^2-y^2)-2xy^2]/(x^2+y^2)^2 whose contourmap of hexagonal symmetry is shown in the accompanying graph. We can interpret this result as representing an Re<<1 viscous 2D flow produced by 6 doublets spaced at 60 degree intervals about the origin. As another case, consider low Re number Couette flow between concentric rotating cylinders. The solution development is shown HERE. I also give you HERE a detailed solution of the biharmonic equation for fluid flow about a small sphere known as the Stokes Drag Formula. Most fluids texts just give the final result but do not go through the detailed steps of the derivation.

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 ftt =c2( fxx+ fyy) , 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*at h=0, with the second equality following from the Bernoulli law for small amplitude waves. We thus find that h oi w =Ak*sinh(kh) and i wA*cosh(kh)=-g h o , from which one obtains at once the well known dispersion relation
w ^2=gk*tanh(kh) which you see plotted in the accompanying graph in non-dimenional form. The corresponding velocity potential of the wave is f =- wh o*cosh[(k(y+h)]*sin( wt-kx)/[k*sinh(kh)] . Note that for shallow water waves, where the wavelength l =2pk is large compared to the water depth h , one finds the very simple and non-dispersive  result  that the velocity equals sqrt[gh] which you are already familiar with from elementary physics. It is a simple matter to show that the fluid elements underneath such waves travels in closed ellliptical orbits so that there will be no net transport of fluid for this linear wave propagation problem.

THE SCHROEDINGER EQUATION AND THE QUANTUM WELL PROBLEM: In  1926, E.Schroedinger converted the classical conservation of energy equation E=p2/(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 Yt, and the potential transfomation V->VY, he arrived at the time-dependent Schroedinger equation i hbar Yt=-(hbar)2/(2m)*(dell2Y)+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. 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(Vo-E))x/hbar]. If one now demands continuity of y and its derivative at x=a, one arrives at the transcendental relation sqrt[(Vo-E)/E]=tan[sqrt(2mEa2/hbar2)]. 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 Vo approaches infinity, the lowest energy level will be at E= p2 hbar2/(8ma2). 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)


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
y(x)=L-1[ (2/s^3)/(1-1/s^2)]=2*[cosh(x)-1]. We have plotted this solution and several of its approximations y[n] obtained by the iteration y[n+1]=x^2+Int[(x-t)*y[n],t=0..x].

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-
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 yn(x) form an orthogonal set with a weight function of r(x) in the range a<x<b. That is Int[r(x)*yn (x)*ym(x), x=a..b]=Const.*d m,n ,where dn,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 yn(x and expand any function f(x) in terms of these as f(x)=Sum[Cn*yn (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 (Fy')] = d/dx[F-y'Fy']=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=n2 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 CYLINDER: -The shortest distance between two points A and B on a surface is referred to as a geodesic. This space curve can be found by extremizing the integral Int[ds]=Int[sqrt(dx^2+dy^2+dz^2)]. For an axisymmetric surface, such as a cone or cylinder, this leads via the Euler-Lagrange equation to the solution theta=a +b*Int[sqrt(1+f'^2)*dr/(r*sqrt[r^2-b^2])], where z=f(r) or r=g(z) and a and b are constants. For the special case of a unit radius cylinder, this geodesic space curve is the helix x=cos(t), y=sin(t), z=t  as shown on the graph.

-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 d4y/dx4 -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>/<y2>, with the bracket notation representing the integral over the indicated range of x. Using the very simple trial function j=(1-x2) 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 yeven=cosh(kx)/cosh(k)-cos(kx)/cos(k) and the odd form yodd=sinh(kx)/sinh(k)-sin(kx)/sin(k), with the eigenvalues l=k4 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)=c1 j1+ c2 j2+c3j3 +--. Here jn represents trial functions each satisfying the imposed homogeneous boundary conditions of the problem. If one now multiplies the Sturm- Liouville equation by jm , substitutes the above series for y(x) , and then integrates over the range a<x<b, there results the equality  Sum[c n (Anm-lB nm ),{n=1..N}]=0, where  Anm=<p(x) j m(x)' jn(x)'+q(x)jn(x) jm(x)> and Bnm=<r(x) jn(x) jm>. 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 cn . For non-trivial solutions for cn, 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 cns by standard algebraic methods involving the minors of the coefficvients.

LOWEST EIGENVALUE AND EIGENFUNCTION FOR THE VIBRATING STRING VIA RAYLEIGH-RITZ:We give here a demonstration of the effectiveness of the Rayleigh-Ritz Method by aplying it to y"+ly=0 subject to y(0)=y(1)=0. Trying the approximation y(x)=c1x(1-x)+c2 x2 (1-x)2, one finds A11 - l B11=1/3-l/30, A12 - lB12=A21 - lB 21=1/15-l/140, and A 22 -B 22=2/105-l/630. So setting the determinant to zero, yields the quadratic equation l 2 /529200-l/4725+1/525=0, which has the two real roots l1=9.86975 and l2=102.13. This result compares very well with the  exact lowest eigenvalue of l1= p2=9.869604 , but not as well with the higher eigenvalue of l 3=9p 2=88.826. The corresponding eigenfunction is found to be y(x)=4x(1-x)[1+0.0127897x(1-x)] which we have plotted in the accompanying graph. Comparing it with the exact lowest eigenfunction of y(x)=sin(px) shows excellent agreement. The above procedure also works for non-selfadjoint differential equations(Galerkin), but there one can no longer be guaranteed that the eigenvalues form an upper bound or will be real.  

APPLICATION OF LAGRANGE MULTIPLIERS: -The Lagrange multiplier method allows one to maximize a function subjected to n constraints. Here one does not deal with integrals but rather with a function W(x,y,z) which is to be extremized subjected to say two constraints f(x,y,z)=0 and g(x,y,z)=0. It is clear that the gradients to these three functions at a fixed spatial point are coplanar, so that grad[W+ l1*f+ l2*g]=0, with  l 1 and l2 being the Lagrange multipliers. This  result leads to three algebaic expressions and these together with the original constraint equations give enough information to solve the problem. Lets demonstrate the procedure by asking -" What is the maximum volume of a rectangular solid which can be placed inside a sphere of radius r=a?" Here the  volume is W=2x*2y*2z and the constraint is x^2+y^2+z^2=a^2. So we have 8yz+ l*2x=0, 8xz+ l*2y, and 8xy+ l*2z=0. These have the obvious solution x=y=z=a/sqrt(3) corresponding to a cube of volume W=8a^3/[3*sqrt(3)]=1.5396 a^3. This maximum rectangular solid volume is only about 37% of the total sphere volume of (4/3)*Pi*a^3. I show you here a Mathematica construction of this cube within a sphere. Putting everything together , including defining the six sides of the cube and partially uncovering the sphere so as to see the cube took a couple hours of effort the first time I attempted this. Try it yourself. It involves using ParametricPlot3D together with the operations Show and Viewpoint, plus some manipulations involving Paintbrush and Graphics Workshop.

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.



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.

    An infinite number of other fractals can also be generated by this program. The particular class of iterations  given by Z[n+1]=Z[n]^m+c yields  figures with m-fold symmetry whose shape changes with the values of the complex starting number c=a+I*b. To enter the realm of computer art and to see a few of these figures which I have generated  for different values of m and c, click HERE, HERE, HERE, HERE , HERE,and HERE. It  is amazing how such intricate  patterns can be generated via such  a simple iteration procedure.

 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 HERE 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 HEREand 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 HERE, 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 HERE to 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 HEREto 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 HERE 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 HERE 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 HERE 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 HERE 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 these HERE.

SERIES SUMMING 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 HERE.

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 HERE. 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-


                                 [ 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


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 HERE 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 HERE to 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 HERE some discussion of these properties and also introduce some additional Fibonacci like sequences.


GENERATION OF BIHARMONIC FUNCTIONS BY COMPLEX VARIABLE METHODS: It is known that a general solution to the biharmonic equation is Psi=xf(x,y)+yg(x,y)+h(x,y)+j(x,y) where f, g, h and j are harmonic functions satisfying the Laplace equation. Knowing that the real and imaginary parts of a function of a complex varaible are harmonic suggests that one can also use complex variable methods to generate biharmonic functions using complex variable methods. Go HERE to see the development.

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 HERE 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(ct2),t=0..x] which leads to well the known Fresnel integrals C(x) and S(x). We discuss HERE 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 HERE 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 HERE 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 HERE 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 (232582657-1) . We discuss HERE some properties of prime numbers and discuss how one might obtain even larger primes by a known extention of the Mersenne  postulate that 2p -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 HERE a discussion in pdf form of how to generate large prime and large composite numbers.

SPHERICAL TRIGONOMETRY, TERRESTRIAL, AND CELESTIAL TRIANGLES: When discussing geodesics in our section on variational methods above, we briefly showed how to find the distance along a geodesic for a sphere by use of Spherical Trigonometry. Let us extend this discussion  by looking a bit  more at how the basic formulas for sperical triangles on the surface of a sphere are derived and then use these results to calculate the length of great circle distance  between two points on the earth and also discuss how to determine  the sun's location on the celetial sphere at any date and time.  GoHERE to see the details.

GRAPHICAL REPRESENTATION OF EVEN, ODD, AND PRIME NUMBERS: In the above discussion on prime numbers we found that all integers can be represented in a single graph with the even and odd integers falling on separate straight lines. We extend this discussion HERE to show, among other things, that for all Mersenne Primes the ratio (M[n]-7)/8 must be an odd integer. Also go HERE to see additional information on primes.

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 HERE 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 HERE 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!+...=j0(1)=0.84147098.., where j0(x) is the spherical Bessel function of the first kind of order zero. Note that this series also equals sin(1).

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 HERE 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 HERE 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 HERE in closed form some of these integrals. 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 HERE 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 Bn(x) generated by texp(xt)/[exp(t)-1] = sum[(Bn(x) /n!)tn, n=0..infinity] and the Euler polynomialsEn(x) generated by 2exp(xt)/[exp(t)+1]=sum[(En(t)/n!)tn, n=0..infinity]. Go HERE to see our discussion on some of the properties of these functions and corresponding numbers and their applications.

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 HERE 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 HERE 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.


USE OF THE FORMULAS S1=2^n+1, S2=2^n+3, S3=2^n-3 AND S4=2^n-1 TO GENERATE PRIME NUMBERS: It is well known that the formula N=2^n-1 can be used to generate the Mersenne primes for certain values of n. Presently there are 46 of these Merseene primes known and one can expect to see many more in coming years as computing power increases. What is sometimes not clear  is that Mersenne primes form only a very small subset  of primes  which can be generated by more general formulas. As we have already shown in an earlier discussion above(see
GRAPHICAL REPRESENTATION OF EVEN, ODD, AND PRIME NUMBERS) all  odd prime numbers  are given by the formulas P1=8n+1, P2=8n+3, P3=8n-3, and P4=8n-1 for appropriate values of n. A subset of the primes generated by these formulas is given by S1=2^n+1, S2-2^n+3, S3=2^n-3, and S4=2^n-1. Go HERE to see a discussion concerning the values of n for which P and S yield primes and when these remain composite numbers. We also show that 2^(2n+1)+1 will always be a composite number for all positive integer values of n. Thus 2^(123456789)+1 is composite. Good luck on finding its factors other than the obvious one of 3.  Use of supercomputers allowed.

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 HERE 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)2)/dm at m=0.5.

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 HERE 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 HERE 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
, one can represent the infinite series for ln(2) as the sum of the finite series sum[(-1)^(n-1)/n, n=1..N]  plus an integral int[t^(N+1)/(1+t), t=0..1]. Since the value of ln(2) and the integral can be evaluated to high accuracy, one simultaneously obtains a precise value for the finite series.

FINDING LOGARITHMS OF NUMBERS: One can define any number N as bn , where b is the base and n=logb(N) its logarithm. Thus log2(1024)=10, loge(10)=2.302585.., and log10 (36)=1.5563025...  . In the era prior to electronic computers and hand caculators most calculations involving the product and quotients of complicated numerical expressions were carried out by means of  logarithms and extensive tables were constructed and utilized for this purpose. We present here a little about the properties of logarithms  and how one goes about calulating  their values. Go HERE to see the discussion.

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 HERE for these discussions.

PROPERTIES OF THE PRIME NUMBERS  GENERATED BY [2^(2n+1)+1]/3,[2^(2n+1)-3]/5, and [2^(2n+1)+3]/5: In the last few months we have been studying how prime numbers are generated. We have found that there are four basic sets of prime numbers  which are generated by 8n+1, 8n-1, 8n+3 and 8n-3  for particular values of integer n. Restricted subsets of these numbers are the Mersenne primes M(n)=2^(2n+1)-1 which are contained in the odd numbers 8n-1 and the Fermat primes which follow from 8n+1. There are also subsets following from 8n+3, 8n-3 and 8n-1 . These include
W(n)=[2^(2n+1)+1]/3 , K(n)=[2^(2n-3)/5, and L(n)=[2^(2n+1)+3]/5. Go HERE to see some of their properties.

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 HERE modular arithmetic to these integers to obtain 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 HERE 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.We want to look HERE at 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 in=im when n=m mod(4k) and k=0,1,2,3.... and that (1+i)4k   is always real while (1+i)4k+2   is pure imaginary.

THE COMPLEX NUMBER N=[a+ib]^n=[b+ia]^n: An interesting complex number is one where [a+ib] and [b+ia] differ from each other but become equal when each is taken to its nth power. That is N=[a+ib]^n=[b+ia]^n. Here the magnitude of N is [a^2+b^2]^(n/2) so that both versions of the number lie at the same point on a circle of radius |N| within the complex plane. Go HERE to see the conditions on a, b and integer n which defines this subset of the family of all complex numbers.

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 HERE 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 HERE 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=(Arn+B/r (n+1)*Ynm(Theta,Phi) where Ynm(Theta,Phi) are the Spherical Harmonics. We extend HERE the discussion on Y
nmgiven in the above PDE course by detailing some of their additional properties including the ability to use Spherical Harmonics in a double Fourier series expansion for representing any function f(Theta,Phi) on the surface of a sphere.

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 HERE 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.

PRIMES OF THE FORM 8n+1 : In our earlier discussions on prime numbers we have found it convenient to plot them as points representing the intersection of the spiral r=(4/Pi)theta and the diagonal lines y=�x. Here we concentrate on those primes lying along the diagonal in the first quadrant and having the form 8n+1 for certain integer values of n. Go HERE to see a pdf file showing how we go about finding the values of n for which 8n+1 is prime. One of the largest primes of the form 8n+1 which we have been able to find with our PC has 1379 digits.

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 HERE and show, for example, that Sum[1/Phi 2n+1, n=0..infinity] =1 and Phi=2sin(3Pi/10).

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 HERE 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 HERE 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 Pn(x). We develop HERE some of their properties including the derivation of generating formulas, the expansion of functions f(x) in terms of these polynomials, and evaluation of some integrals containing the

SOLUTION OF THE INTEGRAL Int[P2n(x)/(a2+x2), x=0..1] : In a previous note we examined some integrals involving the Legendre polynomials. Among these was one which allows one to quickly find approximations to arctan(1/a) expressible as the ratio of two large rational numbers. We give HERE  further details of how one can use such integrals to quickly find very accurate approximations for Pi.

A METHOD FOR QUICKLY FINDING VALUES OF ARCTAN(1/a) FOR LARGE a: In the previous note we examined integrals of the type Int[P2n(x)/(a2+x2),x=0..1. HERE we continue these discussions by examining approximations for arctan(1/a) obtainable from this Legendre polynomial integral. Among other things one finds that arctan(1/a) is closely given by the quotient 5a(11+21a2)/(3(3+30a2+35a4)) when a>>1.

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 HERE 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 HERE 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 D�RER'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 D�rer's Solid or , as I refer to it, as
D�rer's Cube. We look at some of the properties of this solid HERE including  3D graphs showing both the D�rer 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 HERE to see our dicussion on this polyhedron.

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 HERE 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 HERE 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 HERE 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 HERE 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.


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 HERE 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 HERE  and show, among other things, that ln([1+(1/a)] = [1/(6a)][(8+15a-30a2)/(3-5a2)] for a>>1.

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 HERE 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 HERE 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 HERE 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 demonstrateHERE 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 of the Gaussian  HERE 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.

GRAPHICAL DETERMINATION OF PRIME NUMBERS: We examine a graphical technique based on the use of different period sine waves to determine prime numbers by the elimination of neighboring composite numbers. The procedure involves a regrouping of all odd numbers into sequences starting with the lowest prime numbers P=3, 5, 7, 11, etc and having subsequent terms in each sequence generated by P(3+2k) with k=0,1,2 etc. Go HER
E  for the details. Among other results, we show that the Merseene number N=2^11-1=2047 is composite as well as its neighbors 2041, 2043, 2045, 2049, and 2051. The required number of sine curves needed to guarantee a number N is prime equals about 2sqrt(N)/ln(N). To prove that a number is composite will generally require fewer curves.

INTEGRALS INVOLVING BESSEL FUNCTIONS OF THE FIRST KIND: We have discussed in detail the properties of Bessel Functions earlier. HERE 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 examine HERE 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 HERE. 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 HERE for details of these discussions.

TECHNIQUES FOR EVALUATING CERTAIN DEFINITE INTEGRALS: We show you HERE 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 HERE 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 HERE 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 HERE 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 HERE some of their properties 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 HERE to see the discusion.

SOLUTION OF THE NON-LINEAR ALGEBRAIC EQUATION Y2=a+bX2: We examine the integer and non-integer solutions of the equation y2=a+bx2. After converting the equation into a standard hyperbola, we present a solution in terms of continued fractions and treat three special cases in detail. We then examine those solutions [x,y] where these numbrers are integers. It is shown that under certain conditions the ratio xn+1/xn1 =y1+x1sqrt(b). This fact is used to obtain some large integer solutions of
y2=1+2x2 up to 20 digits long. Go HERE to see the discussion.

EVALUATION OF THE RIEMANN ZETA FUNCTION: We evaluate the Riemann Zeta function Zeta(s)=sum[1/ns, n=1..infinity] for different values of real and complex s. The zeros of Zeta[0.5+i tau]) are discussed and a simple  approximation for Zeta(s) for large real s are presented. An exact evaluation of  sum[ln(n)m/ns, n=1..infinity] in terms of the mth derivative of the Zeta(s) function is also found. Go HERE for details.

RIGHT TRIANGLES AND THE BASE AND PYTHAGOREAN TRIPLETS: We examine the solutions to a2+b2=c2 for all cases where a,b,and c are integers. Eulid's formula for the triplets [a,b,c] is derived and it is shown under which conditions such Pythagrean Triplets are also base triplets. Right triangles inscribed and circumscribed by a circle are also treated and the area coverage fraction found. Go HERE for the discussions.

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 HERE to see some variations on their 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 HERE for the details.

MORE ON PRIME NUMBERS: We discuss some more of the properties of prime numbers. Using our x-y plane numbering scheme, we deal with four distinct groups of primes given by N=8n+1, 8n+3, 8n+5, and 8n+7 for appropriate values of n. Mersenne Primes are shown to be sub-class of the group 7, 23, 31, 47, 71, 79, 103,...lying along the diagonal in the fourth quadrant  and defined by N=8n+7. In binary notation all odd numbers, be they prime or composite, lying along this same diagonal end in 111. Go HERE for the details.

CONSTRUCTING SYMMETRIC FIGURES USING DIFFERENT SUB-ELEMENTS: We show HERE 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. We want HERE to look at some of the more important identities involving these functions. It is shown, for example, that Jn(x)"=(1/4)[Jn+2(x)-2Jn(x)+Jn-2(x)].  Also, a quick way to find the zeros of Jn(x) by an iteration process involving the  ratio of  three Bessel Functions Jn(x), Jn-1(x) , Jn+1(x) is shown.

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 HERE how this works for several different cases.

GEOMETRIC PROPERTIES OF PENTAGRAMS: We derive HERE the various sub-lengths of a standard pentagram. 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 HERE the 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 HERE 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 HERE 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: We look HERE at the Fibonacci Sequence and introduce a variation in which the gn+1 term is given as the sum of N terms  having integer values g1=1, g2=2, g3=3, ...gN=N. The ratios of gn+1/gn  for a given N have different constant values as n gets large. The range for these ratios lies between the Golden Ratio  (1+sqrt(5))/2 for N=2 to a ratio of 2 for N going to infinity. Several identities involving such generalized Fibonacci Sequences are derived.

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 2.18  and it is shown that 1/2=sum[(1/pn) -( 1/pn+1 ),n=1..infinity], where pn is the pth prime. Go HERE for the discussion. It is also shown how any composite number may be broken up into the product of primes using a number tree.

MORE 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 HERE to see the discussion.

SYMMETRIC PATTERNS PRODUCED BY THE ITERATION Z[n+1]=Sum{C2n Z[n]2n, n=0..N} : We examine the symmetric patterns created in the complex Z plane where an iteration involving  a  linear combinations of even powers of Z multiplied by specified constants C  remains bounded as n goes to infinity. Some interesting and quite complex patterns emerge for several different examples . Go HERE for the details.

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 HERE to 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 number components in the 1000 bit range and these are typically obtained by a random number search. We want HERE to consider 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 N1=64269+3^96+4^35+7^89, N2=19683+2^59+3^197+5^241, and N3=312+2^346+3^389+4^231.

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. We want HERE to examine Fermat's Little Theorem and also look at a variation which states that  p is prime if Sum[ kn (a^(p-1)-1), n=2..N] mod p = 0.  We also check out the primality of one of our recently found primes p=153,837+4^235+7^528+11^49. It is 455 digits long.

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 HERE at 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[cnbn^xn ,n=1..k] where the cs,bs,and xs are picked at random to fit the desired digit length and A is then varied until P is found to be prime. Go HERE to see the details of our discussion.

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 HERE to see the details of the procedure and a discussion of its limits. Among several other things, we show that
arctan(1/a) is approximately equal to 5a(11+21a2)/(9+90a2+35a4) .

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 HERE 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. HERE 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 HERE 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)].

A NEW HIGHLY ACCURATE METHOD FOR FINDING APPROXIMATE VALUES OF SIN(a), COS(a), AND TAN(a): We continue our discussions on using an integral technique involving the Legendre Polynomials to find quotient approximations for tan(a), sin(a), and cos(a). After obtaining highly accurate approximations for tan(a) equal to K(n,a), where n is the order of the approximation, we find at n=2 that a good approximation for tan(a)  is  K(2,a)=[a^5-105*a^2+945*a]/[15*a^4-420*a^2+945] in -Pi/2<a<+Pi/2. Once K(n,a) has been determined then we find sin(a) and cos(a) via the approximations K(n,a)/sqrt[1+K(n,a)^2] and 1/sqrt[1+K(n,a)^2], respectively. Go HERE for the discussion.

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 HERE 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 y2=1+Nx2 to determine its integer solution pairs (xn,yn). Both a continued fraction form  and a Taylor series expansion  of this BP equation are given. It is shown how highly acurate approximations of the square root of integers can be obtained by using higher solution pairs of the BP equation. We find for example that  the ratio 665857/470832 yields a nine place accurate approximation to the sqrt(2). Go HERE for the details.

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 HERE for 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 HERE 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 HERE 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 (HERE) deals with arctan approximations and the second (HERE) 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(N0)]p can be expressed as [a sqrt(N)-b], where a and b are two easily determined large integers. Go HERE for the discussion. Among others results, we show how sqrt(2) can be evaluated by the difference of just two terms to yield a 48 digit accurate result when p=32.

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 HERE how it 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.

FACTORING SEMI-PRIMES:  We show HERE how to factor a semi-prime N  into its two components p and q by solving the Diophantine equation  y(x)=[(K+x)2-N]/(K+x), where K is an integer near sqrt(N). The x is varied over the range 0<x<sqrt(N) and the integer pair [x,y] found yields p=K+x and q=y-p. As an example, we find for the semi-prime N=89893 and use of K=300 that [x,y]=[73,132] so that p=373 and q=241.

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/pn)s , n=1..x] converge if s is real and greater than one. In this sum, pn is the nth prime number so that p100=541 and p1000=7919. We also evaluate F(s,x) for complex s and find behaviour similar to that of the standard Zeta Function. The product of all primes up to px are examined and some super-composite numbers , such as 
200560490130=(2) (3) (5) (7) (11) (13) (17) (19) (23) (29) (31), and their powers are generated. Go HERE for the details of the discussion.

POWERS OF REAL AND COMPLEX NUMBERS: We evaluate numbers of the form Np where both N  and p can be complex. Using some modular arithmetic and  a binary power tree we carry out several different calculations. Go HERE for the details. Among other things it is shown that [a+ib)(b+ia)]p will always be a real number when p is even. Numbers such as N=2123 are quickly evaluated by noting that this number is equivalent to 8 x 256 x 655367.


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 HERE 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. We derive HERE the 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 HERE 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 HERE 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-

                                                                                                                                2        2
                                                                                                                            3       4       3
                                                                                                                         4     6        6      4
                                                                                                                      5    8       9       8      5
                                                                                                                   6   10   12      12     10    6

We look HERE at 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 HERE 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 HERE 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)2+f(z)3+...f(z)N where f(z) is any function of  the complex variable z=x+iy. Some interesting results are obtained including the fact that   sum[sin(n*pi/4)/2n, n=0..infinity]=2/(5*sqrt(2)-4) and that 1+3+8+20+48+..+N2N-3=2N-2(N-1). Go HERE for the details.

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 HERE 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 HERE 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.

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)2. One can accelerate this type of search by using the exponent vector notation which has as its elements the powers of the prime numbers used in a factorization of a given number. Thus 40 has vector [3 0 1] and 90 the vector [1 2 1]. Adding the exponents in these two vectors produces V=[4 2 2] with has all even elements and hence is a perfect square . Go HERE for further details on this type of squaring. Among other things we show that the googol=10100 can be described by the vector  V=[100 0 100]. Also we show how the semi-prime  24139 can be factored into 101 and 239 by the V vector method.

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 HERE 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((n1(n2-n1)) mod N, N). Here n1 is any large number starting near the square root of N and n2 is a small integer treated as a running integer until gcd assumes a value different from unity. This  non unit answer will be either p or q. Go HERE for details of the technique  including a demonstration of how the semi-primes N=1232895479 and N=4294967297 can be factored with very little effort.

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 HERE for 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 HERE 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)

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)2nan int[ 1/(an2+x2), x=0..1) with an+1=an+sqrt(an2+1) and a0=1  or a0=sqrt(3). In general an will be non-integer but can be made as large as desired making for a rapidly convergent series for Pi. Go HERE to see the details of the discussion.

MERSENNE NUMBERS AND THE INTEGER SPIRAL: We examine the Mersenne Numbers M[n]=22n+1-1 and show that they all lie along a single diagonal in the fourth quadrant of the z plane when one represents all integers as the intersection points between the radial lines theta=Pi n/4 and the Archimedes Spiral z=n exp i(theta). A small subclass of  the M[n]s are found to be prime including M[30], M[44], M[53], M[63], M[260], M[303], and M[639]. Go HERE for the details of the discussion.

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 HERE 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 Jo(x). Go HERE for the details of the discussion.

MORE ON THE SOLUTION OF THE BRAHMAGUPTA EQUATION: The non-linear Diophantine equation y2=1+Nx2 was first studied by the Indian mathematician Brahmagupta(598-668AD) and shown to have integer solutions [N,x,y]. We examine this equation in more detail HERE and show how one can use various methods to obtain integer triplets. One interesting and rapid solution approach is to set y to a constant and then express y2-1 =Nx2 as powers of primes and read off x from those prime terms whose powers are even.  Thus if y=3480, one finds Nx2=72*592*71. Thus x=413 and N=71.

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 HERE 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. It is our purpose HERE 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)�1 , where n represents certain integer values, will be primes. Also we find an important sub-class of primes given by Q=6 x 2n+1�1 which far exceed both Mersenne and Fermat primes in density. Examples are Q=6 x 22023-1 and Q=6 x 2623 +1. Go HERE for the details of the discussion.

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 HERE 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�1. We examine these primes in more detail HERE. 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.

MORE ON FACTORING LARGE PRIMES: In several earlier articles we have discussed methods for factoring large semi-prime numbers. We want HERE to continue these discussions and show , among other things, that 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.

FACTORING N=pq VIA A DIOPHANTINE EQUATION: We show that factoring the semi-prime N=pq=(2n-1)(2m-1) is equivalent to solving the Diophantine  Equation 2U-V=(N-1)/2, where U=nm and V=n+m. One obtains a factor of N at that point were the number s in the solution V=A+2s produces an integer value for  n=[V+sqrt(V2-4U)]/2. The number A is an input to the problem and can be varied. Go HERE for the details of the discussion.

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 HERE 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 fN for certain values of N.  Mersenne Primes M and Q Primes are also discussed. Go HERE for the details.

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 fN, we have found a simple test to check whether a number is prime or composite. The test is applied to both a Fermat Number and a Mersenne Number. In addition we present a graph in the form of a hexagonal spiral which clearly demonstrates that all prime numbers lie along two diagonal lines emanating from the origin of the spiral. Go HERE for the details.

THE BINOMIAL EXPANSION AND ITS MULTIPLE VARIATIONS: We look at the classic Binomial Expansion (a+b)n=sum((n!\/(k!(n-k)!)a(n-k)bk, k=0..infinity) where a and b are allowed to take on complex values and n need not necessarily be integer. Go HERE for the details. Among other things we obtain the binomial expansion for exp(inx), sqrt(2), and for Si(x).

APPROXIMATIONS FOR FUNCTIONS USING ITERATIVE APPROACHES: We show HERE how one can use iteration methods based on the Taylor Series  to approximate the values of functions such as xp, exp(x), arctan(x), and ln(1+x) at specified values of x. We show that sqrt(2) is given by the iteration x[n+1]=x[n]+(2-x[n]2)/(2x[n]) starting with x[0]=1. Also it is shown that excellent approximations to exp(1) are obtainable using the iteration x[n+1]=x[n]{2-ln(x[n])} subject to x[0]=1. The Golden Ratio is evaluated to 80 places based on just eight iterations for sqrt(5). A result of major potential interest is that arctan(1/N) is quickly evaluated by the iteration x[0]=1/N with x[n+1]=x[n]+(cos(x[n]))2[(1/N)-tan(x[n])]. It allows one to find the value of Pi to 100 place accuracy after just six iterations using the Machin Formula.

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])2-(1/2)sin(2z][n]) subject to z[0]=1/N.
The desired solution will be z[infinity]=arctan(1/N) although excellent approximations will already be found at much lower n. The procedure does require an accurate  knowledge of both sin(x) and cos(x). Go HERE for the details of the iteration method. One of the observations arising from this study is that approximations for Pi using arctan formulas are best accomplished  with single term formulas having  larger values for N. Thus iterating Pi=8arctan(1/(1+sqrt(2)) produces an approximation good to 8o places at z[5].

: We look at the integer solutions to a modified Fermat-Pythagorean Formula an+bm=ck, where a, b, c, n, m, and k are all positive integers. Unlike for the Fermat case where n=m=k and no solutions exist when n>2, we show HERE that there are an infinite number of solutions when at least two of the exponents differ from each other. Simple examples occur for 11+23=32, 24+32=52 and 35+102=73. Expressed in terms of a Fermat-Pythagorean Triplet [an, bm, ck]
we are able to generate an infinite number of solutions starting from a few elementary base triplets. One of these generated triplets reads [4n, 23+2n, (3x2n)2] and works for n=0,1,2,etc. We also show how one can write the Mersenne Prime M=217-1=131071 as M=(768)2-48-1.

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 HERE for the details.


GENERATING TRIGONOMETRIC FUNCTIONS USING N SIDED REGULAR POLYGONS: We show how one can express the values of trigonometric functions in terms of radicals using base triangles constructed from n sided regular polygons. We demonstrate how a regular pentagon generates oblique triangles of sides 1-1-phi, where phi is the golden ratio [1+sqrt(5)]/2.   In addition we show that tan(Pi/12)=[sqrt(3)-1]/[sqrt(3)+1] and sin(Pi/8)=sqrt[sqrt(2)-1]/23/4. Go HERE for the details of the analysis.

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 HERE 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=
13469199089641601294165099996964705019718136187846657 factors into 59421121885698253195157962751 x 226673591177742970257407. Go HERE for the details.

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 =2p-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 HERE for the details.

SEMI-PRIMES AND THE SOLUTION OF Y=SQRT(X2-N): It is shown how the semi-prime N=pq is equivalent to the Diophantine equation y2=x2-N. A graph of this equation is presented and it is then shown how one goes about finding the desired integher solution pair [x1, y1] =[(p+q)/2, (p-q)/2]. Specifically we show how  N=12063943 factors into p=5171 and q=2333. Go HERE 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)=|(p2-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 HERE 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 HERE 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 HERE 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 Duplica
tion 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-1
for n=1,2,3.. where n=1,2,3... Prime numbers of 5 or greater have the property that N mod(6)=1 or 5, however
 this is not sufficient unless the ratio test R=N/(6n
+1) or N/(6n-1) fails to yield an integer value. Any odd number for which N mod(6)=3 will always be a composite. Thus N=456723108745338992349075615353107312371 is a composite. Go HERE for
the details including the procedure for factoring large composite numbers into products of primes.
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 inthe range 1<N<700000
is f(665280) =3.39826. Go HERE for the details including the facts that f(pn)=[1-(1/pn-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(a2-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-
x=208671283132153254989370540184529119739903=348181 x 599318409482864530199438051428794563 is achieved. Go HERE for the details.

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 HERE 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 HERE 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)=ax2+2xy+cy2 +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 HERE 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 HERE 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 HERE 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 HERE 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]=Z0 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 HERE 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 HERE 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-xN. This result was the first time an expression had been found which relates the product of  functions with a finite sum. We present HERE a quick proof of the theorem based on the elementary mathematical concepts known to Cotes and predating both calculus and complex variables.

We combine our earlier results on primes and related topics to show that all primes above N=3 have the form N=6n�1 and must have (6n�1)/p be a non-integer for all primes in 5<p<sqrt(N). In modular arithmetic language any prime must have N mod(6) equal to 1 or 5. A number such as N=2674521 is a composite number since N mod(6)=3. A polar plot shows that all positive integer fall at the intersection of a spiral and six radial lines. All primes lie along the radial lines 6n+1  and 6n+5. Go HERE for the details.