Differential Operators  Harmonic Oscillators
Variational Calculus

Operator Calculus:   Linear Operator Rules

Does Mathematica  apply the Fourier Transform/Derivative Rule Automagically?

FourierTransform[D[f[x],{x,1}],x,k]

It does as of version 5.0! .

Table[FourierTransform[D[f[x],{x,i}],x,k],{i,24}]//MatrixForm

Does  Mathematica apply the rules according to the Fourier Transform being a linear operator?

FourierTransform[a f[x] + b g[x],x,k]

Apparently not--so we make some rules that can be applied. It may be instructive to see how to do this.

Two rules are defined so as to factor out non-dependent variables and apply the distribution rules for linear operators

ConstantRule = FourierTransform[(NoX_.  ) *(  fun_ ), x_ , k_] :> NoX FourierTransform[fun,x,k]/;(FreeQ[NoX,x] && MemberQ[fun,x,Infinity])

DistributeRule=  
FourierTransform[Plus[expr_],x_, k_] :>Distribute[FourierTransform[expr,x,k],Plus]

Try them out:

FourierTransform[a g[x]  f[x] ,x, k]//.ConstantRule

FourierTransform[a x f[x] + b v[x]  g[x] + d p[x],x,k]//.DistributeRule//.ConstantRule

Fourier Transforming the Damped and Forced Harmonic Oscillator

Here is the second-order ODE for a damped harmonic oscillator

ODE2nd = Mass D[y[t],{t,2}] + Viscosity D[y[t],t] + SpringK y [t]

Let's  Fourier Transform the left-hand  side of a second-order ODE:

FrrODE2nd =Factor[
FourierTransform[ODE2nd,t,ω]//.DistributeRule//.ConstantRule
]

And now Fourier Transform the right-hand side for a prototype "forced" oscillator with an arbitary frequency ω0. (i.e., K y''  + η y' + m y = cos(ω0 t)

rhs = FourierTransform[Cos[ ω0  t] ,t,ω]

Using the Fourier Transform to Solve the Harmonic Oscillator and Back-Transforming

All the derivatives have been removed, and we can solve for the Fourier transformed as a function of the transformation variable (in this case, ω).

ftsol = Solve[FrrODE2nd ==rhs,FourierTransform[y[t],t,ω]]

Build in some physical assumptions for mass-spring-dashpot systems.

DampedHOAssumptions=  ω0 > 0 &&  Mass > 0 && Viscosity > 0 && SpringK > 0;

Back-transform to find the solution.

FullSimplify[InverseFourierTransform[FourierTransform[y[t],t,ω]/.Flatten[ftsol],ω,t],Assumptions -> DampedHMAssumptions]

Otherwise, the complete solution (i.e., the inhomogeneous plus the homogeneous parts) can be calculated by the DSolve function

GenSol = DSolve[Mass D[y[t],{t,2}] + Viscosity D[y[t],t] + SpringK y[t] == Cos[(Subscript[ω, o]) t],y[t],t]

FullSimplify[y[t]/.Flatten[GenSol],Assumptions->DampedHOAssumptions]

Functionals: An Introduction to Variational Calculus by way of Variation of Parameters: Approximating the Geodesic

The introductory form of calculus of variations involves finding an unknown function, y = y(x),  which maximizes or minimizes a specified definite integral of y(x).

As an introductory problem, we might as a question such as, "of all paths that traverse a rolling landscape, which one is the shortest?"  The y(x) that minimizes the total length is called the geodesic.

The total length of a path over a fixed surface can be written as  an integral of the unknown y(x) (this is worked out in the accompanying lecture notes)
L[y(x)] = (∫_path)ds= (∫_(x=0)^(x=1))(Sqrt[1 + ((((dy/dx)))^2) +((((dh/dx)) )^2)]) dx
Suppose that the "hill" has the shape h(x) = (x^2).

Here, the "real" problem of finding the geodesic from the infinite set of all paths connecting y(x=0) = 0 to y(x=1) = 1  will be replaced by looking at all second-order polynomials: y(x) = a + bx + c(x^2) and then minimizing over the (constrained) parameters. This will be compared with the exact calculation which we will work out below via calculus of variations.

We specify a landscape and an arbitrary cubic path (a function of three parameters)

h =x^2;
$Assumptions= b ∈ Reals;
YGen = a + b x + c (x^2);

Illustrate this surface and identify the starting and ending points.

TheSurface=Show[Plot3D[h,{x,0,1},{y,0,1},PlotStyle-> Opacity[0.5]],Graphics3D[{PointSize[0.05],Darker[Green],Point[{{0,0,0}}],Red,Point[{{1,1,1}}]}]]

The general path must satisfy the boundary conditions:

YBCs =YGen/.(Solve[{(YGen/.x->0)==0,(YGen/.x->1)==1},{a,c}]//Flatten)

There is one remaining free variable, it can be determined by minimizing the integral

TotalDistanceCubic = FullSimplify[Integrate[Sqrt[1 + (D[YBCs,x])^2 + (D[h,x])^2],{x,0,1}]]

To get an idea of how the path' s distance compares with its shape, we can manipulate the path and display that path' s length

Interactive Path Modification and Length Calculation

Manipulate[thedistance= TotalDistanceCubic/.b->B;Show[TheSurface,ParametricPlot3D[{x,YBCs/.b->B,h},{x,0,1},PlotStyle->Thick],Graphics3D[{Text["Distance=" <> ToString[thedistance],{1/2,1/2,1}]}],BaseStyle->Large],{{B,.5},-2,8,.01}]

Minimizing Parameters for the Geodesic's Cubic Approximation

Does the total length of the cubic approximation have a minimum with respect to its remaining parameter?

Plot[TotalDistanceCubic,{b,-2,6}]

OK, the total distance is a minimum somewhere on the interval 0 < b < 1.  Solving for the value of b at the minimum:

BminsolGeodesicCubic= FindMinimum[TotalDistanceCubic,{b,0,1}]

Now we can use the minimizing b to find the best cubic approximation to the geodesic.

GeodesicCubicSolution = YBCs/.BminsolGeodesicCubic[[2]]

This is the shape of the path, it doesn't differ very much from a straight line.

GeodesicCubicPlot = Plot[GeodesicCubicSolution,{x,0,1},PlotStyle-> Thick]

Comparing the Approximation to the Exact Solution

The exact minimizing path can be found by using Calculus of Variations (demonstrated below). The solution is obtained from a boundary-value problem which we do not take up at this point, but it is interesting to see the exact solution and compare it with the approximate one we obtained above.  The closed-form expression for the function that minimizes the climbing time is:

GeodesicExact=((2 x (Sqrt[1+4 (x^2)])+ArcSinh[2 x])/(2 (Sqrt[5])+ArcSinh[2]))

Graphical Comparisons

Plot the exact solution and compare it to our calculated approximation...

GeodesicExactPlot= Plot[GeodesicExact,{x,0,1},PlotStyle->{Thickness[0.005],Hue[1]}];

GraphicsRow[{Show[GeodesicCubicPlot,GeodesicExactPlot] , Show[ Plot[GeodesicExact-GeodesicCubicSolution,{x,0,1},PlotStyle->{Thickness[0.005],Hue[1]}]]}]

These two curves are very close, comparinfg three distances: the cubic approximation, the exact solution, and a straight line.

Distance[f_] := Integrate[Sqrt[1 + (D[f,x])^2 + (D[h,x])^2],{x,0,1}]

Distance[GeodesicExact]

Distance[GeodesicExact]<Distance[GeodesicCubicSolution] < Distance[x]

Approximating the Brachiostone

The "quickest" path of running up-hill  can be estimated with a model of how fast one runs up hill.  Given a starting point y(x=0) = 0 and an ending point y(x=1) = 1, the question is to find, among all possible choices of path y(x) that connects the end and starting point, which is the one that takes the least time if there is a hill h(x) that has to be "run-up."
A reasonable model is that running speed is proportional to Cos(climbing-angle)--that is the velocity is 1 for running on a flat surface and goes to zero as the climbing-angle goes to (π/2) .

Assuming that v(s) = Cos α(s) where α is the angle of incline (i.e., cos α= ((((dx^2)+ (dy^2)))^(1/2))/((((dx^2) + (dy^2) +  (dh^2)) )^(1/2)) then the time increment is
derived from ds/dt= v(s): dt = ds/v(s), then
 time= ∫(ds/(v(s)))=  ∫ ((((dx^2) + (dy^2) +  (dh^2)) )/(((((dx^2)+ (dy^2)))^(1/2)) ))=  ∫ (((1 + ((((dy/dx)))^2) +  ((((dh/dx)))^2)) )/((((1+ ((((dy/dx)))^2)))^(1/2)) ))dx

TotalTimeCubic = FullSimplify[Integrate[(1 + D[YBCs,x]^2 + D[h,x]^2)/Sqrt[1 + D[YBCs,x]^2],{x,0,1}],Assumptions-> b≠1]

Visualizing the Approximation to the Brachiostone

TheSurfaceWithGeodisic= Show[TheSurface,ParametricPlot3D[{x,GeodesicExact,x^2},{x,0,1},PlotStyle->{Thick,Dashed}]]

As with the geodesic, illustrate the dependence of the approximating path for the brachistone and the total time.

Manipulate[thetime = TotalTimeCubic/.b->B;Show[TheSurfaceWithGeodisic,ParametricPlot3D[{x,YBCs/.b->B,x^2},{x,0,1},PlotStyle->{Thick,Red}],Graphics3D[{Text["Time=" <> ToString[thetime],{1/2,1/2,1}]}],BaseStyle->Large],{{B,.5},-2,8,.01}]

Plot[TotalTimeCubic,{b,-2,2}]

BminsolBrachioCubic= FindMinimum[TotalTimeCubic,{b,0,1}]

The minimizer for b can be used to produce the cubic approximation to the brachiostone and then plotted.

BrachioCubicSolution = YBCs/.BminsolBrachioCubic[[2]]

BrachioCubicPlot = Plot[BrachioCubicSolution,{x,0,1},PlotStyle-> Thick]

Calculus of Variations: Euler Equation for the Geodesic

As derived in the notes, the condition for a curve to minimize an integral is given by the variational derivative, setting the variational derivative equal to zero produces the Euler equation

Needs["VariationalMethods`"]

DistanceIntegrand =Sqrt[(1 + (D[y[x],x])^2 + (D[h,x])^2)]

Calling VariationalD from the package "VariationalMethods"

VariationalD[DistanceIntegrand,y[x],x]

Setting the variational derivative to zero is the Euler equation

DistanceExtremalCondition =EulerEquations[DistanceIntegrand,y[x],x]

The Euler equation is an ODE; this one has a closed-form solution.

DistanceMinimizingFunction =DSolve[{DistanceExtremalCondition, y[0]==0,y[1]==1},y[x],x]

DistanceYExactSolution=y[x]/.DistanceMinimizingFunction[[1]]

Calculus of Variations: Numerical Solution of th Euler Equation for the Brachiostone

This is the form of the total time's integrand as derived above for v(s) = cos(α)

TimeIntegrand =(1 + (D[y[x],x])^2 + (D[h,x])^2)/Sqrt[1 + (D[y[x],x])^2]

TimeExtremalCondition =EulerEquations[TimeIntegrand,y[x],x]

This ODE doesn't have a closed-form solution; so we find a numerical approximation to the solution to the Euler equation, extract it and then plot it.

BrachioMinimizerNumerical= NDSolve[{TimeExtremalCondition,y[0]==0,y[1]==1},y[x],x]

BrachioNumerical= y[x]/.BrachioMinimizerNumerical[[1]]

Comparing the Brachiostone to its Approximations

BrachioExactPlot=Plot[BrachioNumerical,{x,0,1},PlotStyle-> {Thick,Darker[Green]}]

GraphicsRow[{Show[BrachioCubicPlot,BrachioExactPlot] , Show[ Plot[BrachioNumerical-BrachioCubicSolution,{x,0,1},PlotStyle->{Thickness[0.005],Hue[1]}]]}]

Time[f_]:=Integrate[TimeIntegrand/.y'[x]-> D[f,x],{x,0,1.0}]

N[Time[BrachioNumerical]]<Chop[Time[BrachioCubicSolution]]<Time[x]