1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
|
# ex_pendulum.ppl
#
# The code in this file is part of Pyxplot
# <http://www.pyxplot.org.uk>
#
# Copyright (C) 2006-2012 Dominic Ford <coders@pyxplot.org.uk>
# 2008-2012 Ross Church
#
# $Id: ex_pendulum.ppl 1261 2012-07-11 21:38:05Z dcf21 $
#
# Pyxplot is free software; you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation; either version 2 of the License, or (at your option) any later
# version.
#
# You should have received a copy of the GNU General Public License along with
# Pyxplot; if not, write to the Free Software Foundation, Inc., 51 Franklin
# Street, Fifth Floor, Boston, MA 02110-1301, USA
# ----------------------------------------------------------------------------
# This script produces a figure for the Pyxplot Users' Guide
# Initialise
reset
# Set output destination
title = "ex_pendulum" ; load "examples/fig_init.ppl"
# BEGIN
set unit angle nodimensionless
theta_approx(a,t) = a*sin(2*pi*t)
theta_exact (a,t) = 2*asin(sin(a/2)*jacobi_sn(2*pi*t,sin(a/2)))
set unit of angle degrees
set key below
set xlabel r'Time / $\sqrt{g/l}$'
set ylabel r'$\theta$'
omega = unit(30*deg)
plot [0:4] theta_approx(omega,x) title 'Approximate solution', \
theta_exact (omega,x) title 'Exact solution'
# END
# Set output destination
load "examples/fig_end.ppl"
title = "ex_pendulum2" ; load "examples/fig_init.ppl"
# BEGIN
subroutine pendulumDivergenceTime(omega, deviation)
{
for t=0 to 20 step 0.05
{
approx = theta_approx(omega,t)
exact = theta_exact (omega,t)
if (abs(approx-exact)>deviation) { ;break; }
}
return t
}
set key top right
set xlabel r'Amplitude of swing'
set ylabel r'Time / $\sqrt{g/l}$ taken to diverge'
set samples 40
plot [unit(5*deg):unit(30*deg)][0:19] \
pendulumDivergenceTime(x,unit(20*deg)) title r"$20^\circ$ deviation", \
pendulumDivergenceTime(x,unit(10*deg)) title r"$10^\circ$ deviation", \
pendulumDivergenceTime(x,unit( 5*deg)) title r"$ 5^\circ$ deviation"
# END
# Call common cleanup script
load "examples/fig_end.ppl"
|
