Higher-Order Ordinary Differential Equations

Forward Differencing Equations for Second-Order Iteration

Recall how we used forward finite differencing to simulate the behavior of a function that changed proportional to its current size in Lecture 19.  Consider now a function that changes "how fast it changes" proportional to its current change and its current value--in other words its current acceleration is proportional to its velocity and to its size:

Let's work out how to find a finite differencing function:

This is the current change or approximation to velocity

ChangePerΔ[F_, i_, Δ_] := (F[i+1] - F[i])/Δ

This is the current change in change or approximation to acceleration

ChangeofChangeperΔ[F_, i_, Δ_] := 
Simplify[(ChangePerΔ[F,i+1,Δ]/Δ - ChangePerΔ[F,i,Δ])]

Our model is that the acceleration is proportional to  size of the current function and its velocity, let these proportions be: -α and -β

DifferenceRelation = ChangeofChangeperΔ[F,i,Δ] == -β  ChangePerΔ[F,i,Δ]- α F[i]

Solve this for the "latest" values

ForDiffSol =Collect[ Solve[DifferenceRelation,F[i+2]],Δ]

Replace to find the form of the solution

ForDiffSolV2 = ForDiffSol/.i-> j-2

Iteration Explicit Sequences of Second-Order Equations

We write a function that takes two values (the current, and its antecedent) and uses a discrete step Δ to push the solution into the future. This method uses the current value to estimate the right-hand-side; so it is an Explicit method.

GrowList[ValuesList_List, Δ_ , α_ , β_] :=Module[{Minus1=ValuesList[[-1,2]], Minus2 = ValuesList[[-2,2]], LastX = ValuesList[[-1,1]]},
Append[ValuesList,{LastX + Δ,
2*Minus1 - Minus2 + Δ *(β*(Minus2 - Minus1) -α*Δ*Minus2)}]]

Notice that this forward differnce needs TWO consectutive past values before it can calculate the current value, create a function to Appends the current value to a list by using the last two values:

The fact that TWO consecutive values are needed is not surprising. Because this iteration procedure is giving an approximation to the solution of a second-order ODE, we need to make a numerical approximation to a second derivative. THREE points—the two consecutive ones and the current one—are needed to do this. (If only two points were used, you could only compute the slope of the straight line connecting the points; if three points are used, a circle can be fit to the three points and the inverse of its radius gives the curvature.)

InitVals = {{0,1},{.001,1}};

result = GrowList[InitVals,.001,1,.1]

Generate a sequence of length 20 from initial values {1,1} for Δ=0.001, α=1, β=0.1

result = GrowList[result,.001,1,.1]

Using NestWhile, we can grow the list until the accumulation of the discrete steps reaches a particular value.

NestWhile[GrowList[#,.001,1,.1]&,InitVals,(Last[#][[1]] < 0.02)&]

Visualizing Second-Order Iteration

Visualize results out to x=50 (takes about 24 seconds to visualize )

ListPlot[NestWhile[GrowList[#,.001,2,.1]&,InitVals,(Last[#][[1]] < 50)&]]

Change parameters for Growth Function (this example shows that the numerical solution does not converge to the accurate solution):

ListPlot[NestWhile[GrowList[#,0.01,2,0]&,InitVals,(Last[#][[1]] <50)&]]

"Textbook" Solutions to Second-Order Homogenous ODEs with Constant Coefficients

We write a general form for the left-hand-side of the homogeneous ODE with constant coefficients β and γ.

TheODE[function_, var_] := D[function[var],{var,2}] + β D[function[var],var] + γ function[var]

TheODE[y,x]

Guess a solution and substitute it into the left-hand side of the ODE:

TheGuess[x_] := Exp[λ x]

TheODE[TheGuess,x]

This will be a solution when the resulting quadratic expression in λ is equal to 0:

λSolution = Solve[TheODE[TheGuess,x]==0,λ]

Set the two solutions to λMinus and λPlus

{λMinus,λPlus}= λ/.λSolution

The linear combination of the two solutiions is the general solution, here we set the two arbitrary constants.

GenSol[x_] := C[LPlus] Exp[λPlus x] +  C[LMinus] Exp[λMinus x]

TheODE[GenSol,z]

This will demonstate that our general solution is, in fact, a solution.

Simplify[TheODE[GenSol,z]]

Building a Visualization of Solution Behavior

Assuming that β and γ are real, the solutions can be characterized by the roots of the quadratic equation. That is, our solution is known to be (e^(λt)):  the solution monotonically grows or shrinks in time if the roots are real-positive or real-negative. The condition that the roots are real is γ <(( (β^2))/4)If the roots have an imaginary part, the solutions can oscillate. Whether the oscillation amplitudes grow or shrink depends on whether the real part of λ is positive or negative. We will calculate conditions that will indicate the nature of solutions:  if the solution is complex, the imaginary part will determine the frequency of the oscillatory solutions. In all cases, the sign of the real part will indicate whether the amplitude will grow exponentially (and become unbounded) or shrink asymptotically towards zero.

RealsCond=Reduce[λPlus∈ Reals &&  λMinus ∈ Reals,{β,γ},Reals]

The Complex and Real Domains

We use RealsCond and RegionPlot to indicate domain where solutions are real or complex.

CplexReal = RegionPlot[RealsCond,{β,-1,1},{γ,-0.25,0.5},AxesLabel->{"β","γ"},PlotStyle-> {Lighter[Blue],Opacity[0.5]}]

Annotate the regions and update the assignment to a new combined graphic.

CplexReal=Show[CplexReal,Graphics[Text["Complex\nCongugate\nRoots",{0,0.4}]],Graphics[Text["Real Roots",{0.0,-0.2}]]]

The Sign of the Real Part of Complex Solutions

If the roots complex, then the sign of the real part is determined only by the sign of β; i.e., λ = ((-β)/2)± (Sqrt[((((β/2)))^2)- γ])

CplexPosNeg =ParametricPlot[{0,t},{t,0,.5},PlotStyle->{Thickness[0.015],Hue[0],Opacity[0.3]}]

CplexPosNeg=Show[CplexPosNeg,Graphics[Text[StyleForm["Positive\nReal\nPart",FontColor->Hue[0]],{-.5,0.3}]],Graphics[Text[StyleForm["Negative\nReal\nPart",FontColor->Hue[0]],{.5,0.3}]]]

CplexPlot =Show[CplexReal,CplexPosNeg]

Real Solutions, Both Positive or Both Negative

Find the conditions where the real roots are positive and/or negative

RealRootsPosCond=Reduce[{λPlus >0,λMinus>0}]

RealRootsPos=RegionPlot[RealRootsPosCond,{β,-1,1},{γ,-0.25,0.5},AxesLabel->{"β","γ"},PlotStyle-> {Darker[Yellow],Opacity[0.5]}]

RealRootsPos= Show[RealRootsPos,Graphics[Text[StyleForm["Positive Roots",FontColor->Hue[.6]],{-0.7,0.025}]]]

AnnotePosRealRoots=Graphics[Text[Style["Positive Roots",FontColor->Hue[0.6]],{-1.,0.025},{-1,0}]];

RealRootsNegCond=Reduce[{λPlus <0,λMinus<0}]

RealRootsNeg=RegionPlot[RealRootsNegCond,{β,-1,1},{γ,-0.25,0.5},AxesLabel->{"β","γ"},PlotStyle-> {Darker[Orange],Opacity[0.5]}]

RealRootsNeg=Show[RealRootsNeg,Graphics[Text[StyleForm["Negative Roots",FontColor->Hue[.6]],{0.7,0.025}]]]

Real Solutions, Opposite Signs

MixedRealRootsCond=Reduce[{λPlus >0,λMinus<0}]

MixedRealRoots=RegionPlot[MixedRealRootsCond,{β,-1,1},{γ,-0.25,0.5},AxesLabel->{"β","γ"},PlotStyle-> {Darker[Blue],Opacity[0.5]}]

MixedRealRoots=Show[MixedRealRoots,Graphics[Text[StyleForm["Mixed Real Roots",FontColor->RGBColor[0,0,.5]],{0.0,-0.05}]]]

Put it all together for a complete graphic....

Show[CplexPlot,RealRootsPos,MixedRealRoots,RealRootsNeg,ImageSize-> Large]

Solving for Boundary Conditions for Second-Order Homogenous ODEs with Constant Coefficients

Second order ODEs require that two conditions  to generate a particular solution.

GenSol[x]

Example of y(0) = 0 and y(L)=1:

SolutionOne = Solve[{GenSol[0]==0,GenSol[L]==1},{C[LPlus],C[LMinus]}]

Write the resulting solution:

SpecificSolutionOne = Simplify[GenSol[x]/.SolutionOne]

Second example:  y(0) = 1 and y'(0)=0

DGen = D[GenSol[x],x]

Now the constants CPlus and CMinus are found by solving:

SolutionTwo = Solve[{GenSol[0]==1,(DGen/.x->0)==0},{C[LPlus],C[LMinus]}]

SpecificSolutionTwo = Simplify[GenSol[x]/.SolutionTwo]

Higher-Order ODEs: the Beam Equation

Set up the ODE solution (we will specify a unit beam stiffness) EI = 1

BeamEquation[y_ ,x_ , w_, BC1_ , BC2_] := DSolve[{y''''[x] == w[x], BC1, BC2},y[x],x]

Functions for Typical Boundary Conditions and Beam Loading

Set up some typical boundary conditions:

A clamp fixes the position and the angle of a beam

Clamp[y_ , x_ , position_, slope_] := {y[x] == position, y'[x] == slope}

The position of a beam and the applied moment (or torque) is applied with a Knob

Knob[y_, x_, position_, moment_] := {y[x]== position, y''[x] == moment}

FreeBeam[y_ , x_] := {y''[x]==0, y'''[x]==0}

Set up some typical beam loads

noload[x_] := 0

A unit load is a force of 1/length applied over the beam length

unitload[x_] := 1

midload appliec a fixed load at x=1/2.

midload[x_] := -10DiracDelta[x - 1/2]

Boxload simulaties a "box" of a fixed width and weight per unit width placed at x

boxload[x_,pos_,width_,mag_] := Piecewise[{{mag,pos-width/2 ≤ x ≤ pos+width/2}}]

Simulated linearly decreasing distributed load (pushes up when x<1/2; down for x > 1/2)

linearload[x_] :=500*(1/2 - x)

BeamEquation[y,x,noload,Clamp[y,0,0,0],Knob[y,1,-.1,0]]

Visualizing Solutions to the Beam Equation

Plot[Evaluate[
y[x]/.BeamEquation[y,x,noload,Clamp[y,0,0,0],Knob[y,1,-.25,0]]
],{x,0,1}]

BeamViz: Visualization Function

BeamViz[DistLoadx_, BC1_, BC2_] :=  
Plot[Evaluate[
y[x]/.BeamEquation[y,x,DistLoadx,BC1,BC2]]
,{x,0,1},PlotRange->{-0.5,0.5},AspectRatio->1,PlotStyle->{Thickness[0.03],Hue[0]}]

Unit load pushing up, a clamp at the left and a free left end: This is a first-order model for the deflection of an airplane wing.

BeamViz[unitload,Clamp[y,0,0,0],FreeBeam[y,1]]

Fixed Knobs that rotate freely at both ends, load: Approximates the defllection of a sieve loaded by  shearing flow.

BeamViz[linearload,Knob[y,0,0,0],Knob[y,1,0,0]]

A box of width 1/8 sitting at x=3/4 for a double cantilever beam.

BeamViz[boxload[#,3/4,1/8,-500]&,Clamp[y,0,0,0],Clamp[y,1,0,0]]

A point load at the center of a beam clamped at the left and sitting on a pinned point at the right.

BeamViz[midload,Clamp[y,0,0,0],Knob[y,1,0,0]]

A load applied to the center with balancing loads applied at the ends: Three Point Loading

BeamViz[midload,Knob[y,0,0,0],Knob[y,1,0,0]]

Frivolous Visualization Example: the Diving Board Beam

We wish to simulate the deflection of a diving board as a diver walks toward the end.  A diving board may be modeled as beam with constant cross-section. The boundary conditions are that the board is clamped at the beginning; has a pivot located somewhere near the center; and the "diving" end . For "competitive" diving boards, the pivot point is adjustable. Here we fix the pivot at x=1/2.

This next function is a silly little graphic for the walking diver.

Graphics Functions for Animating Board and Diver

humanbean[theta_] :=Module[{redtheta=N[Pi/8 Sin[theta]]},
{Pink,
Rotate[Disk[{0,0.05},{0.005,0.05}],redtheta,{0,0.1}],
Rotate[Disk[{0,0.125},{0.004,0.035}],2redtheta,{0,0.16}],
Darker[Pink],Disk[{0,0.12},{.015,.05}],Disk[{0,0.17},0.015], 
Rotate[Disk[{0,0.05},{0.005,0.05}],-redtheta,{0,0.1}],
Rotate[Disk[{0,0.125},{0.004,0.035}],-2redtheta,{0,0.16}],

Darker[Red],Rectangle[{-0.015,0.06},{0.015,0.09}],
Darker[Red],Text["MiT",{-0.009,.075},BaseStyle->{FontSize-> 3,FontFamily-> "Helvetica"},Background-> Gray]
}];

Manipulate[diver=humanbean[t];Graphics[Translate[diver,{1,1}]],{t,0,10}]

The general solution for this beam equation is a bit complicated.  We approximate the diver's force with a Delta-function. To get a solution, we have to match boundary conditions at the center-pivot.  If the diver is to the left of the pivot, we use the implied slope at the pivot as input to the slope at the pivot for the beam to the right of the pivot.  If the diver is to the right of the pivot, we allow the slope to be an undetermined parameter, and then enforce that the shear force is continous across the pivot. The result is a fairly involved piece-wise function.

board[loc_] := board[loc] =Module[{yleft,midslope,yright,y,mund},
yleft=yright=0;
If[0 ≤  loc ≤ 0.4999,
yleft=y[x]/.Flatten[DSolve[{y''''[x]==-8DiracDelta[x-loc],y[0]==0, y'[0]==0, y[1/2]==0,y''[1/2]==0},y[x],x]];
midslope = D[yleft,x]/.x->1/2;
yright= y[x]/.Flatten[DSolve[{y''''[x]==0,y''[1]==0, y'''[1]==0, y[1/2]==0,y'[1/2]==midslope},y[x],x]]
];
If[0.5001≤  loc ≤ 1,
yright= y[x]/.Flatten[DSolve[{y''''[x]==-8DiracDelta[x-loc],y''[1]==0, y'''[1]==0, y[1/2]==0,y'[1/2]==mund},y[x],x]];
midslope = D[yright,x]/.x->1/2;
yleft=y[x]/.Flatten[DSolve[{y''''[x]==0,y[0]==0, y'[0]==0, y[1/2]==0,y'[1/2]==mund},y[x],x]];
midslope = mund/.Flatten[Solve[(D[yleft,{x,2}]/.x->1/2)== (D[yright,{x,2}]/.x->1/2),mund]];
yleft= yleft/.{mund-> midslope};
yright = yright/.{mund->midslope}];
Piecewise[{{yleft,0 < x <1/2},{yright,1/2<x ≤ 1}}]
]

This is a fixed graphic for the pivot.

pivot = Graphics[{Darker[Green], Polygon[{{1/2,0},{1/2-.02,-0.1},{1/2+ .02,-0.1}}], Black,Disk[{1/2,0},.01],Gray,Disk[{1/2,0},.005]}];

And finally, our interactive simulation for a the walking diver.

Manipulate[
boardsol = N[ board[Rationalize[pos,.001]]];
diver =humanbean[10 π pos];
plot = Plot[boardsol,{x,0,1},PlotRange-> {-.6,0.2},ImageSize->Large,PlotStyle->{Thickness[0.025],Darker[Blue]},AspectRatio->8/10];
ypos=N[ (boardsol/.x->pos)]/.HeavisideTheta[0.]->1/2;
Column[{Show[plot,Graphics[Translate[diver,{pos,ypos}]],pivot],boardsol}],
{{pos,2/7,"position"},0,1-.001,.001}]