Simulering av inverterad pendel med analytisk mekanik

INOM

EXAMENSARBETE TEKNIK, GRUNDNIVÅ, 15 HP

STOCKHOLM SVERIGE 2019 ,

Simulated inverted

pendulum using analytical mechanics

CARL CLAUSON

MALTE SLOMA

INOM

EXAMENSARBETE TEKNIK, GRUNDNIVÅ, 15 HP

STOCKHOLM SVERIGE 2019 ,

Simulering av inverterad pendel med analytisk

mekanik

CARL CLAUSON

MALTE SLOMA

Contents

1 Abstract 2

2 Sammanfattning 3

3 Introduction 4

3.1 Background . . . . 4

3.2 Purpose . . . . 4

3.3 Problems . . . . 5

3.3.1 Double Pendulum . . . . 5

3.3.2 Kapitza Pendulum . . . . 6

3.3.3 Kapitza Pendulum with elliptic oscillation. . . . . 7

4 Theory 8 4.1 Newtonian Mechanics . . . . 8

4.2 Functionals and Variational Calculus . . . . 8

4.3 Lagrange’s equations . . . . 10

5 Method 11 5.1 Software . . . . 11

5.1.1 Maple with Sophia package . . . . 11

5.1.2 Matlab . . . . 11

5.2 General method . . . . 11

5.3 Double pendulum . . . . 12

5.3.1 Lagrangian mechanics . . . . 12

5.3.2 Solving with Euler-Lagrange equations . . . . 13

5.4 Kapitza pendulum . . . . 14

5.4.1 Lagrangian mechanics . . . . 14

5.4.2 Solving with Euler-Lagrange equations . . . . 14

5.5 Kapitza pendulum with elliptical oscillation . . . . 15

5.5.1 Lagrangian mechanics . . . . 15

5.5.2 Solving with Euler-Lagrange equations . . . . 15

6 Results 16 6.1 Double pendulum . . . . 16

6.2 Kapitza pendulum . . . . 17

6.2.1 Minimum angular velocity to maintain stability . . . . 17

6.2.2 Maximum initial angle to maintain stability . . . . 20

6.3 Kapitza pendulum with elliptical oscillation . . . . 22

6.3.1 Maximum width of the orbit to maintain stability . . . . 22

6.3.2 Influence of width in macro oscillations . . . . 24

A Appendix code 29 A.1 Double pendulum . . . . 29

A.2 Kapitza’s pendulum . . . . 30

A.3 Kapitza’s pendulum with elliptical oscillation . . . . 31

1 Abstract

Analytical mechanics is an alternative to Classical (Newtonian) mechanics for

calculating the movement of particles and systems of particles. In this paper we

analyze three different mechanical systems using analytical mechanics. We then

implement this into the maple based software ”Sophia” to further analyze how

changing some of the parameters alter the behavior of the system. One of our

main focuses in the report is studying the stability of an inverted pendulum. In

1951 the Russian physicist Pyotr Kapitza, figured out that the pendulum could

be stabilized in the upside-down position by forcing it into a rapid vertical oscil-

lation. We examine the boundaries of this stability and explore if the pendulum

can instead be stabilized using an elliptical oscillation.

2 Sammanfattning

Analytisk mekanik ¨ ar ett alternativ till klassisk (Newtonsk) mekanik f¨ or att

ber¨ akna r¨ orelsen av partiklar och system av partiklar. I den h¨ ar rapporten anal-

yserar vi tre olika mekaniska system med analytisk mekanik. Vi implementerar

sedan detta i den maple-baserade programvaran ”Sophia” f¨ or att ytterligare

analysera hur ¨ andringar av vissa parametrar f¨ or¨ andrar systemets beteende. Ett

av v˚ ara huvudsyften med den h¨ ar rapporten ¨ ar att studera stabiliteten hos en

inverterad pendel. ˚ Ar 1951 fann den ryska fysikern Pyotr Kapitza att pendeln

kunde stabiliseras i dess upp-och-ner-l¨ age genom att tvinga den till en snabb

vertikal sv¨ angning. Vi studerar gr¨ anserna f¨ or denna stabilitet och unders¨ oker

om pendeln ist¨ allet kan stabiliseras genom en elliptisk sv¨ angning.

3 Introduction

3.1 Background

Newtonian mechanics has throughout its existence been used to determine vec- tor quantities of forces, momenta and motion in a certain system. When scien- tists discovered that the Newtonian mechanics was not enough to solve certain types of mechanical problems, analytical mechanics was invented which can be defined as a collection of closely related alternative formulations of Newtonian mechanics. Analytical mechanics utilizes scalar properties of motion, i.e usually the systems total potential energy and kinetic energy and not vector quantities of accelerations and forces of individual particles. Lagrangian mechanics is an analytic method which is ideal for calculating the motion of different systems with various types forces, for example pendulums. In the present work we will restrict the forces acting on the system to conservative fields.

A regular pendulum is in it’s ground state stable and will all the time strive to keep the center of mass as low down as possible. In 1951 the Russian physicist Pyotr Kapitza, figured out that the pendulum could be stabilized in the upside- down position by forcing the pendulum into a quick vertical oscillation. This pendulum is one of threee systems that will be studied in this thesis.

3.2 Purpose

The methods of analytical mechanics will be applied to model the three different

physical systems using the dynamic tool for Maple; Sophia [2]. To compute the

motion of these systems using regular Newtonian mechanics is possible but even

with only two degrees of freedom it becomes very complicated. Therefore will

the selected systems firstly be modeled quite easily using analytical mechanics

and then use these methods to solve and analyze three different problems. The

systems are a regular double pendulum, and then two different versions of the

Kapitza pendulum [3].

3.3 Problems

3.3.1 Double Pendulum

Problem 3.29 from Apazidis [3].

Figure 1: Double pendulum.

’Consider a double pendulum consisting of two bars OA and AB of masses m

and m

and lengths l

and l

respectively that are free to rotate in the vertical plane according to the figure 1. Choose the following numerical values of the parameters m

= 3 kg, m

= 1 kg, l

= l

= 1 m and calculate and plot the trajectory of end B by means of the Sophia and Graphics packages. Choose the following initial conditions q

(0) = 1 rad, ˙ q

(0) = 0 rad, q

(0) = 1.5 rad,

˙

q2(0) = 0 rad. Calculate then the trajectory of the point B for slightly different

initial conditions and show the sensitive dependence of the motion of the system

on initial conditions by comparing the two trajectories.’

3.3.2 Kapitza Pendulum

Figure 2: Kapitza pendulum.

Consider an inverted pendulum consisting of one bar OA of mass m

and length l

that is free to rotate in the vertical plane according to the figure. The bar will be forced into a harmonic vertical oscillation which will stabilize the pendulum.

In order to change the behaviour of the pendulum, parameters will be changed and examined.

The stability of the pendulum will be studied by deciding the minimum angular

velocity ω and the maximum initial angle to maintain stability. The parameters

of the pendulum is described in the method section.

3.3.3 Kapitza Pendulum with elliptic oscillation.

Figure 3: Kapitza pendulum with elliptic oscillation.

Consider an inverted pendulum consisting of one bar OA of mass m

and length l

that is free to rotate in the vertical plane according to the figure. The bar will be forced into a harmonic oscillation around an ellipse which will stabilize it vertically. In order to change the behaviour of the pendulum, the parameters will be changed and examined.

The stability of the pendulum will be studied by deciding the maximum re-

lationship between the height and the width of the ellipse in order to maintain

vertical equilibrium. The effect of an increasing width will also be studied,

causing macro oscillations.

4 Theory

4.1 Newtonian Mechanics

From Newton’s laws of motions, the concept of energy can be derived. When describing a conservative mechanical system, energy becomes a key factor to solve the often non-linear differential equations for the system. In classical mechanics the energy of the system can be divided into two parts. The kinetic energy is defined as

T = m(v · v)

2 (4.1)

for a particle with all of it’s mass distributed in one single point. When talking about larger and more complicated structures than particles, e.g a bar, the kinetic energy will consist of two parts, translational and rotational:

T

= m(v

· v

)

2 + ω

Iω

2 (4.2)

Where the inertia I and the angular velocity ω per definition Apazidis [3].

Potential energy, is defined per the Cambridge Dictionary as the ”energy stored by something because of position (as when an object is raised), because of its condition (as when something is pushed or pulled out of shape), or in chemical form (as in fuel or an electric battery)” [5].

4.2 Functionals and Variational Calculus

A functional is a mathematical operator that uses functions as its input argu- ments and returns a scalar. One fundamental problem of variational calculus is to find a real function y(x) of a real variable x so that the functional of I[y]

becomes an extreme value.

I[y] = Z

x1

f y(x), y

(x), x
dx (4.3)

I[y] is the functional of y with fixed endpoints x

and x

. The purpose of the functional is to determine those functions y(x), which uses the given values y

= y(x

) and y

= y(x

) as endpoints and makes the functional I[y] an extremum. In other words, for which functions y(x) the functional I[y] assumes a maximum, minimum or a saddle point. Now assume

I[β] = Z

x1

f y(x, β), y

(x, β), x
dx (4.4)

where y(x, β) = y(x) + βη(x) with η(x

) = η(x

) = 0. This implies that y(x) is different from y(x, β) but have the same boundary conditions, i.e. multiple paths are possible between the two endpoints.

Figure 4: Example of functions with same boundary conditions.

The variation of I is given by

δI = dI dβ dβ =

Z

x1

dx ( ∂f

∂y dy dβ + ∂f

∂y

dy

dβ

)

dβ. (4.5)

When the second term in the integral is integrated by parts Z

x1

dx ∂f

∂y

d dx

dy dβ

!

= − Z

x1

dx dy dβ

d dx

∂f

∂y

! + ∂f

∂y

dy dβ     

x2

x1

(4.6)

the conclusion is that the boundary terms do not contribute because dy/dβ = η(x). Therefore

δI = Z

x1

dx ( ∂f

∂y − d dx

∂f

∂y

) dy

dβ dβ. (4.7)

The expression

∂f

∂y − d dx

∂f

∂y

= δf

δy (4.8)

is called the variational derivative of f by y. δy can be seen as a small variation of the curve y(x) since (dy/dβ)dβ = δy. When δI = 0 i.e. I(β) assumes a extreme value, the integrand in (4.7) can be negligible. This results in

∂f

∂y − d dx

∂f

∂y

= 0 (4.9)

which is Euler’s differential equation of variational calculus.

4.3 Lagrange’s equations

The Lagrange function is defined as the difference of the systems kinetic and potential energy.

L = T − V (4.10)

As implemented, the Lagrange function is a function of generalized coordinates q

, their time derivatives ˙ q

and the time t

L = L(q

, ˙ q

, t) (4.11)

where the index k varies depending on the degree of freedom. The action integral S is defined as an integral of the Lagrange function (4.11) along a path in the configuration space between two fixed times t

and t

.

S = Z

t1

L(q

, ˙ q

, t) (4.12)

The action integral (4.12) can be used in Hamilton’s Variational Principle which is a central concept for Langranian mechanics. The principle states that the path of the system between the times t

and t

has a extremum along the actual path.

This implies that the action intregal calculated using the actual path is less than for any another possible path, if the extremum is a minimum. Applied to all possible paths, the Hamilton’s Variational Principle can be rewritten using variation calculus from Section 2

δS = δ Z

t₁

L(q

, ˙ q

, t) = 0 (4.13)

and with this derive the Euler-Lagrange Equation which is the equation used in Sophia and throughout this thesis.

∂L

∂q

− d dt

∂L

∂ ˙ q

!

= 0 (4.14)

5 Method

5.1 Software

To solve the non-linear differential equations described, different kinds of soft- ware were used to solve and present the given problems. The softwares used will be described briefly and the areas of application.

5.1.1 Maple with Sophia package

Maple is a computer software system, that provides symbolic solutions of math- ematical problems and technical calculations [2]. Sophia is a package custom made for Maple, which both simplifies the construction of a mechanical system and the method of understanding the system.

The main perks of the Sophia package is that time dependent variables can be declared. If the variables are declared in different coordinate systems, a relationship between the coordinate systems can be established. Sophia also provides a function which makes it easy to evaluate different vectors in whichever already declared coordinate system. These functions are crucial in order to construct the description of the mechanical systems for this thesis (see Appendix A for each system).

5.1.2 Matlab

Matlab is a computer software mainly used for visualization of functions through graphs and numerical computing. To present the data from the different mecha- nial systems, graphs and animations were created by importing data from Maple to Matlab.

5.2 General method

The three different mechanical systems in this thesis were analyzed using the same procedure. This model describes a general solution for each system. The specific equations for each system is described later in this section.

I. By examining the mechanical system, the generalized coordinates, con- stants and relations between coordinate systems are defined.

II. Assemble the Lagrangian L = T − V , where T is the kinetic energy and V the potential energy of the system.

III. Determine and calculate the Euler-Lagrange equations.

IV. Animate and plot the generalized coordinates and their respective deriva-

tives numerically using Maple’s ode-solver.

5.3 Double pendulum

The two bars OA and AB have masses m

and m

and lengths l

and l

re- spectively. Since there is no friction in the joins O and A the pendulum can move freely in the xy-plane. The system can be considered conservative due to that the only acting force on the pendulum is gravitation. The generalized coordinates q

and q

are defined as the angle between the two bars OA and AB respectively.

Figure 5: Double pendulum.

5.3.1 Lagrangian mechanics

The Lagrangian L is defined as the difference between the kinetic and the po- tential energy of the system, from equation (4.10). The kinetic energy of the bar OA is defined as

T

= m

L

ω

3 (5.1)

where ω

is the angular velocity. Since the bar OA is fixed in O, the contri-

bution to the kinetic energy of the system will only depend on the rotational

motion. With the parallel axis theorem, equation (5.1) can be derived. How-

ever, the bar AB will experience both translative and rotational motion. The

kinetic energy for the bar AB is

T

= m

L

ω

+ ω

6 + m

(u

· u

)

2 (5.2)

where u

is the speed of the bar AB. The total kinetic energy is then calculated

T

= T

+ T

= m

L

ω

3 + m

L

ω

+ ω

6 + m

(u

· u

) 2 . (5.3)

The reference level of the system is the x-axis of the N -frame and the potential energy of the system is therefore calculated by the vertical displacement of the centre of masses with regard to this level. The potential energy of the bar OA is

V

= − m

gL

cos(q

)

2 (5.4)

and for the bar AB

V

= gm

(−L

cos(q

) − 2L

) cos(q

)

2 + L

sin(q

) sin(q

) 2

!

. (5.5)

The total potential energy is then calculated

V

= − m

gL

cos(q

)

2 +gm

(−L

cos(q

) − 2L

) cos(q

)

2 + L

sin(q

) sin(q

) 2

! . (5.6) With equation (4.10) the Lagrangian L can be calculated using the defined energies in this section.

5.3.2 Solving with Euler-Lagrange equations

Using the theory from the previous section, the Euler-Lagrange equations can be

derived. The ode-function of Maple is then used to numerically solve the system

for given initial conditions. Afterwards, plots can be designed and presented.

5.4 Kapitza pendulum

The bar OA have mass m

length l

. Since there is no friction in the joint O, the pendulum can move freely in the xy-plane. The system can be considered conservative due to that the only acting force on the pendulum is gravitation.

The generalized coordinate q

is defined as the angle between the y-axis and the pendulum in an upright position.

Figure 6: Kapitza pendulum. The bar OA oscillates vertically.

5.4.1 Lagrangian mechanics

The Lagrangian is defined as the difference between the kinetic and the potential energy of the system, from equation (2.8). The kinetic energy of the bar is defined as

T = m(u

· u

)

2 + m

L

ω

3 (5.7)

The reference level of the system is the x-axis of the N-frame and the potential energy of the system is therefore calculated by the vertical displacement of the centre of mass with regard to this level. The potential energy of the bar

V = m

g A sin(ωt) + cos(q

)L

2

!

(5.8)

5.4.2 Solving with Euler-Lagrange equations

Using the theory from the previous section, the Euler-Lagrange equations can be

derived. The ode-function of Maple is then used to numerically solve the system

for given initial conditions. Afterwards, plots can be designed and presented.

5.5 Kapitza pendulum with elliptical oscillation

The bar OA have mass m

length l

. Since there is no friction in the joint O, the pendulum can move freely in the xy-plane. The system can be considered conservative due to that the only acting force on the pendulum is gravitation.

The generalized coordinate q

is defined as the angle between the y-axis and the pendulum in an upright position.

Figure 7: Kapitza pendulum with elliptic oscillation. The bar OA oscillates elliptically.

5.5.1 Lagrangian mechanics

The Lagrangian is defined as the difference between the kinetic and the potential energy of the system, from equation (2.8). The kinetic energy of the bar is defined as

T = m(u

· u

)

2 + m

L

ω

3 (5.9)

The reference level of the system is the x-axis of the N-frame and the potential energy of the system is therefore calculated by the vertical displacement of the centre of mass with regard to this level. The potential energy of the bar

V = m

g A sin(ωt) + cos(q

)L

2

!

(5.10)

5.5.2 Solving with Euler-Lagrange equations

Using the theory from the previous section, the Euler-Lagrange equations can be

derived. The ode-function of Maple is then used to numerically solve the system

for given initial conditions. Afterwards, plots can be designed and presented.

6 Results

6.1 Double pendulum

(a) Initial condition q

= 1,0 q

= 1,4. (b) Initial condition q

= 1,1 q

= 1,5.

Figure 8: Trajectory plots for the pendulum from given initial conditions for the first 10 seconds of the simulation.

(a) Initial condition q

= 1,0 q

= 1,4. (b) Initial condition q

= 1,1 q

= 1,5.

Figure 9: Phase portraits for the pendulum from given initial conditions for the

first 100 seconds of the simulation.

6.2 Kapitza pendulum

6.2.1 Minimum angular velocity to maintain stability

The first study of this pendulum is to find the smallest angular velocity of the pendulum ω to maintain stability. The initial angle is set to q

= 10

rad since smaller initial angles have a smaller effect on the pendulums trajectory. The graphs in this section describes the angle q

(y-axis) over time (x-axis).

Figure 10: ω = 36 rad/s. Not stable.

The pendulum is close to upright to begin with but turns unstable quickly. Since

the initial angle is very small, one can see that small changes could have big

impact for larger t.

Figure 11: ω = 36, 3 rad/s. Not stable but falls down slower.

When we increase the frequency to 36, 3 rad/s, the pendulum falls over much

slower compared to the frequency of 36 rad/s.

Figure 12: ω = 36, 324 rad/s. Stable pendulum.

The pendulum has now reached the required oscillation frequency required for

stability. This means that the critical ω is less than 37 rad/s and bigger than

36 rad/s. Our closest approximation of 36, 324 rad/s can be seen in figure 12.

Figure 13: ω = 37 rad/s. Stable and increased frequency of maximum. The pendulum is stable for all ω > 36, 324 rad/s.

Discussion:

In this study we set out to find the minimum angular velocity required to keep the pendulum stable. To mimic reality where a small disturbance would cause the pendulum to swing slightly from its state of equilibrium we initiated the simulation with the pendulum at a very slight angle. Specifically q

= 10

. By gradually increasing the angular velocity of the oscillation we can see in Figure (10), (11) and (12) the pendulum is unstable all the way up to a frequency of ω = 36, 324, after which it is stable.

6.2.2 Maximum initial angle to maintain stability

The second study of this pendulum is to calculate the largest initial angle pos-

sible in order to maintain stability.

Figure 14: q

= 0, 5 rad. Stable pendulum.

Figure 15: q

= 0, 521 rad. Stable pendulum.

Figure 16: q

= 0, 522 rad. Not stable pendulum.

Discussion:

In this study we wanted to find what the maximum initial angle while maintain- ing stability is. The frequency used for this was very high at ω = 1000 rad/s.

This was to make sure that any instability was not caused by a lack of angular velocity, By gradually increasing the initial angle of the oscillation we can see in Figure (14), (15) and (16) the pendulum is stable all the way up to an initial angle of 0, 521 rad, which is approximately 29, 9 °, after which it is unstable.

1,55 rad, approximately 88.8 °, is the maximum angle of the pendulum which can be seen in figure (15).

6.3 Kapitza pendulum with elliptical oscillation

6.3.1 Maximum width of the orbit to maintain stability

The first study of this pendulum is to decide the maximum width of the orbit,

with a given fixed height A = 0, 1.

Figure 17: B = 0, 098 m. Stable pendulum.

Figure 18: B = 0, 09854 m. Stable pendulum.

Figure 19: B = 0, 09855 m. Not stable pendulum.

From these graphs we can see that when using a very high frequency ω = 1000 rad/s and very small initial angle q

= 10

m the pendulum is stable for oscillation widths of 0, 09854 m and lower.

Discussion:

When using an elliptical oscillation pattern instead of a vertical one we noticed that the pendulum is unstable for big enough widths of the ellipse. To find the breaking point between stability and instability we increased the width gradu- ally until the system was no longer stable. A high frequency (ω = 1000 rad/s) and small initial angle (q

= 10

) was used to minimize their affect on the maximum possible width. We found that the pendulum is stable for oscillation widths of 0, 09854 m and lower.

6.3.2 Influence of width in macro oscillations

Here we are considering a Kapitza pendulum with the frequency of ω = 50 rad/s

and an initial angle of 0, 2 rad.

Figure 20: B = 0 m.

Figure 21: B = 0, 03 m.

Figure 22: B = 0, 05 m.

Discussion:

When studying the graphs of different Kaptiza pendulums it is clear that while

the base is oscillating quickly with a low amplitude, the pendulum is oscillating

slowly and with a larger amplitude around it’s equilibrium point. In this study

we looked at the impact on these macro oscillations of changing the width of

the base oscillation, specifically with the frequency of ω = 50 rad/s and an

initial angle of 0, 2 rad. From Figure 20-22 we can see that by making the base

oscillation more elliptical the frequency of the macro oscillations decreases.

References

[1] Kojevnikov, A. (2019). Pyotr Leonidovich Kapitsa. [online] Encyclopæ- dia Britannica. Available at: https://www.britannica.com/biography/Pyotr- Kapitsa [Accessed 30 Apr. 2019]

[2] Maple. (2019). Ontario, Canada. University of Waterloo.

[3] Apazidis, N. (2012). Mekanik 2, Partikelsystem, stel kropp och analytisk mekanik. Lund: Studentlitteratur.

[4] Apazidis, N. (2012). Mekanik 1, Statik och partikeldynamik. Lund: Stu- dentlitteratur.

[5] Cambridge Dictionary. Potential energy. [online] Available at:

https://dictionary.cambridge.org/dictionary/english/potential-energy [Ac- cessed 02 May 2019].

[6] Scheck F. Mechanics : From Newton’s Laws to Deterministic Chaos. Berlin,

Heidelberg : Springer Berlin Heidelberg; 2010.

A Appendix code

A.1 Double pendulum

> restart;

> with(plots): with(plottools):

> read ‘C:/Users/calle/OneDrive/Skrivbord/Sophia/SophiaV6N.txt‘;

> read ‘C:/Users/calle/OneDrive/Skrivbord/Sophia/Graphics.txt‘;

> dependsTime(q1, u1, q2, u2);

> &rot [N, A, 3, q1];

> &rot [A, B, 3, q2];

> r1 := A &ev [0, -L1/2, 0];

> r11 := (2 &** r1);

> r2 := r11 &++ (B &ev [0, -L2/2, 0]);

> r22 := r11 &++ (2 &**(B &ev [0, -L2/2, 0]));

> v1 := N &fdt r1;

> v2 := N &fdt r2;

> v1 := simplify(v1);

> v2 := simplify(v2);

> T := 1/3*m1*L1^2*(q1t^2) + 1/2*(1/12*m2*L2^2+m2*L2^2/4)*(q1t+q2t)^2+1/2*m2*(v2 &o v2);

> simplify(T);

> V := m1*g*(N &to r1)[1][2] + m2*g*(N &to r2)[1][2];

> L := T - V;

> L := subs(q1t=u1, q2t=u2, L);

> L_q1 := diff(L, q1);

> L_u1 := diff(L, u1);

> L_u1_t := &dt L_u1;

> eq1 := L_u1_t - L_q1 = 0;

> eq1 := subs(q1t=u1, q2t=u2, eq1);

> L_q2 := diff(L, q2);

> L_u2 := diff(L, u2);

> L_u2_t := &dt L_u2;

> eq2 := L_u2_t - L_q2 = 0;

> eq2 := subs(q2t=u2, q2t=u2, eq2);

> kde := {q1t=u1, q2t=u2};

> eq1 := {u1t = solve(eq1,u1t)};

> eq2 := {u2t = solve(eq2,u2t)};

> eqs := eq1 union eq2 union kde;

> InitCond := {q1(0)=1.1, u1(0)=0, q2(0)=1.5, u2(0)=0};

> eqst := subs(toTimeFunction, eqs);

> param := {m1=3 ,m2=1, L1=1, L2=1, g=9.8}: eqst := subs(param, eqst);

> ff := dsolve(eqst union InitCond, {q1(t), q2(t), u1(t), u2(t)}, type=numeric, maxfun=500000);

> odeplot(ff,[q2(t),u2(t),color=blue],0..100,thickness=1,numpoints=1000);

> #animatepath_2d(ff,r22,param,0,10,100,100,red,2,N);

> plotpath_2d(ff,r22,param,0,10,10000,blue,2,N)

A.2 Kapitza’s pendulum

> restart;

> with(plots): with(plottools):

> read ‘C:/Users/calle/OneDrive/Skrivbord/Sophia/SophiaV6N.txt‘;

> read ‘C:/Users/calle/OneDrive/Skrivbord/Sophia/Graphics.txt‘;

> dependsTime(q, u);

> &rot [N, A, 3, q];

> r1 := N &ev [0, b*sin(omega*t), 0];

> r2 := r1 &++ (A &ev [0, L1/2, 0]);

> v1 := simplify(N &fdt r1);

> v2 := simplify(N &fdt r2);

> T := 1/2*m*(v2 &o v2) + 1/24*m*L1^2*qt^2;

> V := m*g*(N &to r2)[1][2];

> L := T - V;

> L := subs(qt=u, L);

> L_q := diff(L, q);

> L_u := diff(L, u);

> L_u_t := &dt L_u;

> eq := L_u_t - L_q = 0;

> eq := subs(qt=u, qt=u, eq);

> kde := {qt=u};

> eq := {ut = solve(eq,ut)};

> eqs := eq union kde;

> InitCond := { q(0)=0.2, u(0)=0};

> eqst := subs(toTimeFunction, eqs);

> param := {m=1, L1=1, g=9.8, b=0.1, omega = 50}: eqst := subs(param, eqst);

> ff := dsolve(eqst union InitCond, {q(t), u(t)}, type=numeric, maxfun=500000);

> odeplot(ff,[[t,q(t)]],0..5,thickness=0,numpoints=1000);

> #animatepath_2d(ff,r2,param,0,20,1000,200,red,2,N);

A.3 Kapitza’s pendulum with elliptical oscillation

> restart;

> restart;

> with(plots): with(plottools):

> read ‘C:/Users/calle/OneDrive/Skrivbord/Sophia/SophiaV6N.txt‘;

> read ‘C:/Users/calle/OneDrive/Skrivbord/Sophia/Graphics.txt‘;

> dependsTime(q2, u2);

> &rot [N, A, 3, q2];

> r1 := N &ev [a*cos(omega*t), b*sin(omega*t), 0];

> r2 := r1 &++ (A &ev [0, L1/2, 0]);

> v1 := simplify(N &fdt r1);

> v2 := simplify(N &fdt r2);

> T := 1/2*m2*(v2 &o v2) + 1/24*m2*L1^2*q2t^2;

> V := m2*g*(N &to r2)[1][2];

> L := T - V;

> L := subs(q2t=u2, L);

> L_q2 := diff(L, q2);

> L_u2 := diff(L, u2);

> L_u2_t := &dt L_u2;

> eq2 := L_u2_t - L_q2 = 0;

> eq2 := subs(q2t=u2, q2t=u2, eq2);

> kde := {q2t=u2};

> eq2 := {u2t = solve(eq2,u2t)};

> eqs := eq2 union kde;

> InitCond := { q2(0)=0.2, u2(0)=0};

> eqst := subs(toTimeFunction, eqs);

> param := {m2=1, L1=1, g=9.8, a=0.05, b=0.1, omega = 50}: eqst := subs(param, eqst);

> ff := dsolve(eqst union InitCond, {q2(t), u2(t)}, type=numeric, maxfun=500000);

> odeplot(ff,[[t,q2(t)]],0..5,thickness=0,numpoints=1000);

> #animatepath_2d(ff,r2,param,0,1,1000,200,red,2,N);