Question 1.10.4: Compare the responses of the nonlinear and linear pendulum e......
Compare the responses of the nonlinear and linear pendulum equations using numerical integration and the value (g / \ell)=(0.1)^2, for (a) the initial conditions x_0=0.1\ rad and v_0=0.1\ rad / s, and (b) the initial conditions x_0=1\ rad and v_0=1\ rad / s, by plotting the responses. Here x and v are used to denote \theta and its derivative, respectively, in order to accommodate notation available in computer codes.
Depending on which program is used to integrate the solution numerically, the equations must first be put into first-order form, and then either Euler integration or Runge–Kutta routine may be implemented and the solutions plotted. Integrations in Matlab, Mathematica, and Mathcad are presented. More details can be found in Appendix F. Note that the response to the linear system is fairly close to that of the full nonlinear system in case (a) with slightly different frequency, while case (b) with larger initial conditions is drastically different. The Mathcad solution follows.
First, enter the initial conditions for each response:
| v _0:=1 \quad x _0:=1 \quad v 1_0:=0.1 \quad x 1_0:=0.1 |
Next, define the frequency and the number of and size of the time steps:
| \omega:=0.1 \quad N :=2000 \quad i :=0 . . N \quad \Delta:=\frac{6 \cdot \pi}{\omega \cdot N } |
The nonlinear Euler integration is
| \left[\begin{array}{l} x_{i+1} \\ v_{i+1} \end{array}\right]:=\left[\begin{array}{c} v_i \cdot \Delta+x_i \\ -\omega^2 \cdot \sin \left(x_i\right) \cdot \Delta+v_i \end{array}\right] |
The linear Euler integration is
| \left[\begin{array}{l} x 1_{i+1} \\ v 1_{i+1} \end{array}\right]:=\left[\begin{array}{c} v 1_i \cdot \Delta+x 1_i \\ -\omega^2 \cdot(\Delta) \cdot x 1_i 1+v 1_i \end{array}\right] |
The plot of these two solutions yields
(Figure 1)
Here the dashed line is the linear solution. Next, compute these solutions again using initial conditions close to unstable equilibrium values:
| \begin{gathered} x _0:=\pi \quad v _0:=0.1 \quad x 1_0:=\pi \quad v 1_0:=0.1 \\ {\left[\begin{array}{l} x _{ i +1} \\ v _{ i +1} \end{array}\right]:=\left[\begin{array}{c} v _{ i } \cdot \Delta+ x _{ i } \\ -\omega^2 \cdot \sin \left( x _{ i }\right) \cdot \Delta+ v _{ i } \end{array}\right]\left[\begin{array}{l} x 1_{ i +1} \\ v 1_{i+1} \end{array}\right]:=\left[\begin{array}{c} v 1_{ i } \cdot \Delta+ x 1_{ i } \\ -\omega^2 \cdot(\Delta) \cdot x _{ i } 1+ v 1_i \end{array}\right]} \end{gathered} |
(Figure 2)
Script File
The Matlab code for running the solutions (using Runge–Kutta this time) and plotting is obtained by first creating the appropriate M-files (named lin_pend_dot.m and NL_pend_dot.m, defining the linear and nonlinear pendulum equations, respectively).
function xdot = lin_pend_dot(t,x)
omega = 0.1; % define the natural frequency
xdot(1,1) = x(2);
xdot(2,1) = -omega^2*x(1);
function xdot = NL_pend_dot(t,x)
omega = 0.1; % define the natural frequency
xdot(1,1) = x(2);
xdot(2,1) = -omega^2*sin(x(1));
In the command mode type the following:
% Overplot linear & nonlinear simulations of the free
% response of a pendulum.
x0 = 0.1; v0 = 0.1;
ti = 0; tf = 200;
% linear
[time_lin,sol_lin]=ode45('lin_pend_dot',[ti tf],[x0 v0]);
% nonlinear
[time_NL,sol_NL]=ode45('NL_pend_dot',[ti tf],[x0 v0]);
% overplot displacements
figure
plot(time_lin,sol_lin(:,1),'-')
hold
plot(time_NL,sol_NL(:,1),'—')
xlabel('time (s)')
ylabel('theta')
title(['Linear vs. nonlinear pendulum with x0 = ' ...
num2str(x0) ' and v0 = ' num2str(v0)])
legend('linear','nonlinear')
Here the plots have been suppressed as they are similar to those from the Mathcad solution. Next, consider the Mathematica code to solve the same problem.
First, we load the add-on package that will enable us to add a legend to our plot.
In[1]:= <<PlotLegends'
Define the natural circular frequency, ω
In[2]:= ω =0.1;
The following cell solves the linear differential equation, then the nonlinear differential equation, and then produces a plot containing both responses.
In[3] := xlin=NDSolve[{xl''[t]+ω2*xl[t]==0, xl[0]==.1, xl'[0]==.1},
xl[t],{t,0,200}]
xnonlin=NDSolve[{xnl''[t]+ω2*Sin[xnl[t]]==0, xnl[0]==.1, xnl'[0]==.1},
xnl[t], {t,0,200}]
Plot[{Evaluate[xl[t]/.xlin], Evaluate[xnl[t]/.xnonlin]},{t,0,200},
PlotStyle S {Dashing[{}], Dashing[{.01,.01}]}, PlotLabel S "Linear and
Nonlinear Response, Stable Equilibrium",
AxesLabel S {"time,s",""}, PlotLegend S {"Linear","Non-Linear"},
LegendPosition S {1,0}, LegendSize S {.7, .3}]
Out[3]= {{xl[t] S InterpolatingFunction[{{0., 200.}},<>][t]}}
Out[4]= {{xnl[t] S InterpolatingFunction[{{0., 200.}},<>][t]}}
Out[5]=
In[6]:= Clear[xlin, xnonlin, xl, xnl]
xlin=NDSolve[{xl''[t]+ω2*xl[t] == 0, xl[0] ==π, xl'[0] ==.1}, xl[t],
{t, 0, 200}]
xnonlin=NDSolve[{xnl''[t]+ω2*Sin[xnl[t]] ==0, xnl[0] ==π, xnl'[0] ==.1},
xnl[t], {t, 0, 200}]
Plot[{Evaluate[xl[t]/.xlin], Evaluate[xnl[t]/.xnonlin]}, {t, 0, 200},
PlotStyle S {Dashing[{}], Dashing[{.01,.01}]}, PlotRange S {-20, 40},
PlotLabel S "Linear and Nonlinear Response, Unstable Equilibrium",
AxesLabel S {"time,s", ""},
PlotLegend S {"Linear", "Nonlinear"}, LegendPosition S {1, 0},
LegendSize S {.7, .3}]
Out[7]= {{xl[t] S InterpolatingFunction[{{0., 200.}},<>][t]}}
Out[8]= {{xnl[t] S InterpolatingFunction[{{0., 200.}},<>][t]}}
Out[9]=
