### 4.10 Plot the constant energy levels for a nonlinear pendulum

Problem:

Plot the constant energy levels for the nonlinear pendulum in $$\theta ,\dot {\theta }$$

Assume that $$m=1,g=9.8m/s^{2},L=10m$$

The constant energy curves are curves in the Y-X plane where energy is constant. The Y-axis represents $$\dot {\theta }$$, and the X-axis represents $$\theta$$

We assume the pendulum is given an initial force when in the initial position ($$\theta =0^{0}$$) that will cause it to swing anticlock wise. The pendulum will from this point obtain an energy which will remain constant since there is no damping.

The higher the energy the pendulum posses, the larger the angle $$\theta$$ it will swing by will be.

If the energy is large enough to cause the pendulum to reach $$\theta =\pi$$ (the upright position) with non zero angular velocity, then the pendulum will continue to rotate in the same direction and will not swing back and forth.

The expression for the energy $$E$$ for the pendulum is ﬁrst derived as a function of $$\theta ,\dot {\theta }$$\begin {align*} E & =PE+KE\\ & =mgL\left ( 1-\cos \theta \right ) +\frac {1}{2}mL^{2}\dot {\theta }^{2} \end {align*}

This was solved in Mathematica using the ListContourPlot[] command after generating the energy values.

The original Matlab implementation is left below as well as the Maple implementation. However, these were not done using the contour method, which is a much simpler method. These will be updated later.

The following is the resulting plot for $$m=1,g=9.8m/s^{2} ,L=10m$$

Mathematica

SetDirectory[NotebookDirectory[]];
energy[theta_, speed_, m_, g_, len_] := m*g*len*(1 - Cos[theta]) + (1/2)*m*len^2*speed^2;
m = 1; g = 9.8; len = 10;
data = Table[energy[i, j, m, g, len], {j, -4, 4, 0.01}, {i, -3*Pi, 3*Pi, Pi/20.}];
ListContourPlot[data, InterpolationOrder -> 2, Contours -> 20, MaxPlotPoints -> 30,
DataRange -> {{-3*Pi, 3*Pi}, {-4, 4}},
FrameTicks -> {{{-4, -3, -2, 0, 1, 2, 3, 4}, None},
{{-3*Pi, -2*Pi, -Pi, 0, Pi, 2*Pi, 3*Pi}, None}},
FrameLabel -> {{"$Theta]'(t)", None}, {"\[Theta](t)", "constant energy levels for nonlinear pendulum model"}}, ImageSize -> 400, LabelStyle -> 14, RotateLabel -> False]  Matlab function HW2 %HW2 problem. MAE 200A. UCI, Fall 2005 %by Nasser Abbasi %This MATLAB function generate the constant energy level %curves for a nonlinear pendulum with no damping close all; clear; m=1; g=9.8; L=10; %number of data points (positio vs speed) to find per curve nPoints=40; nEnergyLevels = 8; %number of energy levels lowAnglesToVisit = linspace(0,pi,nEnergyLevels); lowEnergyLevels(1:nEnergyLevels)=m*g*L*(1-cos(lowAnglesToVisit)); highAnglesToVisit = linspace(pi,2*pi,2*nEnergyLevels); highEnergyLevels=zeros(1,2*nEnergyLevels); initialHighEnergy=2*m*g*L; Q=0.2; for i=1:2*nEnergyLevels highEnergyLevels(i) = initialHighEnergy+(Q*i*initialHighEnergy); end A = zeros(length(lowAnglesToVisit)+length(highAnglesToVisit),2); A(:,1) = [lowAnglesToVisit highAnglesToVisit]; A(:,2) = [lowEnergyLevels highEnergyLevels]; [nAngles,~]=size(A); data=zeros(nPoints,2); for j=1:nAngles currentAngle=A(j,1); currentEnergyLevel =A(j,2) ; angles=linspace(0,currentAngle,nPoints); data(1:nPoints,1)=angles(:); for m=1:nPoints data(m,2)=speed(currentEnergyLevel,angles(m)); end doPlots(data); end title(sprintf(['constant energy curves, nonlinear pendulum,'... 'quantum=%1.3f'],Q)); xlabel('angle (position)'); ylabel('speed'); set(gca,'xtick',[-4*pi,-3*pi,-2*pi,-pi,0,pi,2*pi,3*pi,4*pi]); set(gca,'XTickLabel',{'-4pi','-3pi','-2pi','-pi',... '0','pi','2pi','3pi','4pi'}) print(gcf, '-dpdf', '-r600', 'images/matlab_e77_1'); function s=speed(energy,angle) m=1; g=9.8; L=10; if angle<pi s=sqrt(2/(m*L^2)*(energy-m*g*L*(1-cos(angle)))); else s=sqrt(2/(m*L^2)*(energy-m*g*L*(1+cos(angle-pi)))); end function doPlots(data) plotCurves(data,0); plotCurves(data,2*pi); plotCurves(data,-2*pi); function plotCurves(data,k) plot(data(:,1)+k,data(:,2)); hold on; plot(data(:,1)+k,-data(:,2)); plot(-data(:,1)+k,-data(:,2)); plot(-data(:,1)+k,data(:,2));  Maple The Maple solution was contributed by Matrin Eisenberg  restart; plotEnergyLevels:= proc(m::positive, g::positive, L::positive, Erange::positive, numContours::posint) local plotOptions, maxTheta, thetaDot, plotContour, Emax, energies, contours; plotOptions := args[6..-1]; maxTheta := E- arccos(max(1-E/m/g/L, -1)); thetaDot := E- theta- sqrt(max(2/m/L^2 * (E - m*g*L*(1-cos(theta))), 0)); plotContour := E- plot(thetaDot(E), 0..maxTheta(E), numpoints=5, plotOptions); # Create first quadrant of middle region of the graph. Emax := Erange*m*g*L; energies := {Emax * n/numContours  n=1..numContours}; if Erange 2 then energies := energies union {2*m*g*L}; fi; contours := map(plotContour, energies); # Create other quadrants. map2(rcurry, plottools[reflect], {[[0,0],[0,1]], [0,0], [[0,0],[1,0]]}); contours := contours union map(f- map(f, contours)[], %); # Create left and right regions of the graph. map2(rcurry, plottools[translate], {-2*Pi,2*Pi}, 0); contours := contours union map(f- map(f, contours)[], %); # Complete the graph. plots[display]([%[], plot([[-2*Pi,0], [0,0], [2*Pi,0]], plotOptions, style=point)], view=[-3*Pi..3*Pi, -thetaDot(Emax)(0)..thetaDot(Emax)(0)], title="Constant energy levels, undamped pendulum", labels=["angle", "angular speed"], xtickmarks=[seq](evalf(k*Pi)=sprintf("%a pi", k), k=-3..3)); end: plotEnergyLevels(1, 9.8, 10, 5, 15); ;  sectionSolve numerically the ODE u””+u=f using point collocation method Problem: Give the ODE \[ \frac {d^{4}u\left ( x\right ) }{dx^{4}}+u\left ( x\right ) =f$ Solve numerically using the point collocation method. Use 5 points and 5 basis functions. Use the Boundary conditions $$u\left ( 0\right ) =u\left ( 1\right ) =u^{\prime \prime }\left ( 0\right ) =u^{\prime \prime }\left ( 1\right ) =0$$

Use the trial function $g\left ( x\right ) ={\displaystyle \sum \limits _{i=1}^{N=5}} a_{i}\left ( x\left ( x-1\right ) \right ) ^{i}$ Use $$f=1$$

The solution is approximated using $$u\left ( x\right ) \approx g\left ( x\right ) ={\displaystyle \sum \limits _{i=1}^{N=5}} a_{i}\left ( x\left ( x-1\right ) \right ) ^{i}$$.

$$N$$ equations are generated for $$N$$ unknowns (the unknowns being the undetermined coeﬃcients of the basis functions). This is done by setting the error to zero at those points. The error being $\frac {d^{4}g\left ( x\right ) }{dx^{4}}+g\left ( x\right ) -f$ Once $$g\left ( x\right )$$ (the trial function is found) the analytical solution is used to compare with the numerical solution.

Mathematica

Clear["Global*"];
nBasis  = 5;
nPoints = 5;
a = Array[v,{nBasis}]


Out[392]= {v[1],v[2],v[3],v[4],v[5]}


trial    = Sum[a[[n]] (x(x-1))^n,{n,1,nBasis}];
residual = D[trial,{x,4}]+trial-1;
mat      = Flatten[Table[residual==0/.x->n/(2*nPoints-1),
{n,1,nPoints-1}]];

mat = Join[mat,{2 a[[1]]+2a[[2]]==0}];
sol = N@First@Solve[mat,a]


Out[397]= {v[1.]->-0.0412493,
v[2.]->0.0412493,
v[3.]->0.000147289,
v[4.]->-0.0000245233,
v[5.]->-3.28259*10^-8}


numericalSolution=Chop@
FullSimplify@Sum[sol[[n,2]](x(x-1))^n,
{n,1,nBasis}]

numericalSolutionDiff4 = D[numericalSolution,{x,4}];
numericalMoment = D[numericalSolution,{x,2}];



\begin {equation*} \begin {split} 0.0412493 x-0.0826459 x^3+0.0416667 x^4 -0.00034374 x^5-1.51283*10^-8 x^6\\ +0.0000984213 x^7-0.0000248515 x^8 +1.64129*10^-7 x^9-3.28259*10^-8 x^10 \end {split} \end {equation*}

(*now analytical solution is obtained
using DSolve and compared to numerical solution*)
eq = u''''[x]+u[x]==1;
bc = {u[0]==0,u[1]==0,u''[0]==0,u''[1]==0};
analyticalSolution=First@DSolve[{eq,bc},u[x],x];
analyticalSolution=u[x]/.analyticalSolution;
analyticalSolution=FullSimplify[analyticalSolution]



$-\frac {\cos \left (\frac {1-(1+i) x}{\sqrt {2}}\right )+\cos \left (\frac {1-(1-i) x}{\sqrt {2}}\right )+\cosh \left (\frac {1-(1+i) x}{\sqrt {2}}\right )+\cosh \left (\frac {1-(1-i) x}{\sqrt {2}}\right )-2 \left (\cos \left (\frac {1}{\sqrt {2}}\right )+\cosh \left (\frac {1}{\sqrt {2}}\right )\right )}{2 \left (\cos \left (\frac {1}{\sqrt {2}}\right )+\cosh \left (\frac {1}{\sqrt {2}}\right )\right )}$

analyticalMoment = D[analyticalSolution,{x,2}];



$-\frac {-i \cos \left (\frac {1-(1+i) x}{\sqrt {2}}\right )+i \cos \left (\frac {1-(1-i) x}{\sqrt {2}}\right )+i \cosh \left (\frac {1-(1+i) x}{\sqrt {2}}\right )-i \cosh \left (\frac {1-(1-i) x}{\sqrt {2}}\right )}{2 \left (\cos \left (\frac {1}{\sqrt {2}}\right )+\cosh \left (\frac {1}{\sqrt {2}}\right )\right )}$

dispError[pt_] :=
((analyticalSolution-numericalSolution)/.x->pt)/
(analyticalSolution/.x->.5)

momentError[pt_]:=
((analyticalMoment-numericalMoment)/.x->pt)/
(analyticalMoment/.x->.5)

(*Now the percentage displacement
and moment errors are plotted*)
Plot[dispError[z],{z,0,1},
ImageSize->400,
Frame->True,
GridLines->Automatic,
GridLinesStyle->Dashed,
PlotStyle->Red,
FrameLabel->{{"error",None},
{"x","percentage displacement error"}},
LabelStyle->14]



Plot[momentError[z],{z,0,1},
ImageSize->400,
Frame->True,
GridLines->Automatic,
GridLinesStyle->Dashed,
PlotStyle->Red,
FrameLabel->{{"error",None},
{"x","percentage moment error"}},
LabelStyle->14,PlotRange->All]



(*Now the Residual error distribution
over x is plotted*)
Plot[Abs[numericalSolutionDiff4+
numericalSolution-1],{x,0,1},
ImageSize->400,
Frame->True,
GridLines->Automatic,
GridLinesStyle->Dashed,
PlotStyle->Red,
FrameLabel->{{"error",None},
{"x","Residual error distribution over x"}},
LabelStyle->14,
PlotRange->All]



Maple

 #Solve u''''+u=1 using point collocation. Zero B.C.
#written by Nasser Abbasi for fun on April 25,2006.
restart;
Digits:=15;
nBasis  := 5:
nPoints := 5:
coef    := [seq(a[i],i=1..nBasis)]:
g := (x,n)-coef[n]*(x*(x-1))^n:  #basis function
f := x-sum(g(x,k),k=1..nBasis):  #trial function
moment := x-diff(f(x),x$2): residual := x-diff(f(x),x$4)+f(x)-1:
A := seq(subs(x=N/(2*nPoints-1),residual(x)=0),N=1..nPoints-1):
A := A,2*(coef[1]+coef[2]=0):
sol   := solve([A],coef):
coef  := map(rhs,sol[1]):
evalf(coef):
numericalSolution:=x-sum(coef[n]*(x*(x-1))^n  ,n=1..nBasis):
The numerical solution is , numericalSolution(x);
numericalSolutionMoment := x-diff(numericalSolution(x),x$2): #Now obtain exact solution from Maple eq := diff(W(x),x$4)+W(x)=1:
bc    := W(0)=0,W(1)=0,D(D(W))(0)=0,D(D(W))(1)=0:

exact := unapply(rhs(dsolve({eq,bc})),x):

exactMoment := x-diff(exact(x),x\$2):

displacmentError := x-(exact(x)-numericalSolution(x))/exact(.5):
momentError := x- (exactMoment(x)-numericalSolutionMoment(x))/evalf(subs(x=.5,exactMoment(x))):

plot(displacmentError(x),x=0..1,labels=['x',"Disp. error %"]);
plot(momentError(x),x=0..1,labels=['x',"Moment error %"]);
`