vcs.maemo.org Git - gnuplot/blob - demo/bivariat.dem

   1 #
   2 # $Id: bivariat.dem,v 1.8 2006/07/06 21:52:19 sfeam Exp $
   3 #
   4 # This demo is very slow and requires unusually large stack size.
   5 # Do not attempt to run this demo under MSDOS.
   6 #
   7
   8 # the function integral_f(x) approximates the integral of f(x) from 0 to x.
   9 # integral2_f(x,y) approximates the integral from x to y.
  10 # define f(x) to be any single variable function
  11 #
  12 # the integral is calculated using Simpson's rule as
  13 #          ( f(x-delta) + 4*f(x-delta/2) + f(x) )*delta/6
  14 # repeated x/delta times (from x down to 0)
  15 #
  16 delta = 0.2
  17 #  delta can be set to 0.025 for non-MSDOS machines
  18 #
  19 # integral_f(x) takes one variable, the upper limit.  0 is the lower limit.
  20 # calculate the integral of function f(t) from 0 to x
  21 # choose a step size no larger than delta such that an integral number of
  22 # steps will cover the range of integration.
  23 integral_f(x) = (x>0)?int1a(x,x/ceil(x/delta)):-int1b(x,-x/ceil(-x/delta))
  24 int1a(x,d) = (x<=d*.1) ? 0 : (int1a(x-d,d)+(f(x-d)+4*f(x-d*.5)+f(x))*d/6.)
  25 int1b(x,d) = (x>=-d*.1) ? 0 : (int1b(x+d,d)+(f(x+d)+4*f(x+d*.5)+f(x))*d/6.)
  26 #
  27 # integral2_f(x,y) takes two variables; x is the lower limit, and y the upper.
  28 # calculate the integral of function f(t) from x to y
  29 integral2_f(x,y) = (x<y)?int2(x,y,(y-x)/ceil((y-x)/delta)): \
  30                         -int2(y,x,(x-y)/ceil((x-y)/delta))
  31 int2(x,y,d) = (x>y-d*.5) ? 0 : (int2(x+d,y,d) + (f(x)+4*f(x+d*.5)+f(x+d))*d/6.)
  32
  33 set autoscale
  34 set title "approximate the integral of functions"
  35 set samples 50
  36 set key bottom right
  37
  38 f(x) = exp(-x**2)
  39
  40 plot [-5:5] f(x) title "f(x)=exp(-x**2)", \
  41   2/sqrt(pi)*integral_f(x) title "erf(x)=2/sqrt(pi)*integral_f(x)", \
  42   erf(x) with points
  43
  44 pause -1 "Hit return to continue"
  45
  46 f(x)=cos(x)
  47
  48 plot [-5:5] f(x) title "f(x)=cos(x)", integral_f(x)
  49
  50 pause -1 "Hit return to continue"
  51
  52 set title "approximate the integral of functions (upper and lower limits)"
  53
  54 f(x)=(x-2)**2-20
  55
  56 plot [-10:10] f(x) title "f(x)=(x-2)**2-20", integral2_f(-5,x)
  57
  58 pause -1 "Hit return to continue"
  59
  60 f(x)=sin(x-1)-.75*sin(2*x-1)+(x**2)/8-5
  61
  62 plot  [-10:10] f(x) title "f(x)=sin(x-1)-0.75*sin(2*x-1)+(x**2)/8-5", integral2_f(x,1)
  63
  64 pause -1 "Hit return to continue"
  65
  66 #
  67 # This definition computes the ackermann. Do not attempt to compute its
  68 # values for non integral values. In addition, do not attempt to compute
  69 # its beyond m = 3, unless you want to wait really long time.
  70
  71 ack(m,n) = (m == 0) ? n + 1 : (n == 0) ? ack(m-1,1) : ack(m-1,ack(m,n-1))
  72
  73 set xrange [0:3]
  74 set yrange [0:3]
  75
  76 set isosamples 4
  77 set samples 4
  78
  79 set title "Plot of the ackermann function"
  80
  81 splot ack(x, y)
  82
  83 pause -1 "Hit return to continue"
  84
  85 set xrange [-5:5]
  86 set yrange [-10:10]
  87 set isosamples 10
  88 set samples 100
  89 set key top right at 4,-3
  90 set title "Min(x,y) and Max(x,y)"
  91
  92 #
  93 min(x,y) = (x < y) ? x : y
  94 max(x,y) = (x > y) ? x : y
  95
  96 plot sin(x), x**2, x**3, max(sin(x), min(x**2, x**3))+0.5
  97
  98 pause -1 "Hit return to continue"
  99
 100 #
 101 # gcd(x,y) finds the greatest common divisor of x and y,
 102 #          using Euclid's algorithm
 103 # as this is defined only for integers, first round to the nearest integer
 104 gcd(x,y) = gcd1(rnd(max(x,y)),rnd(min(x,y)))
 105 gcd1(x,y) = (y == 0) ? x : gcd1(y, x - x/y * y)
 106 rnd(x) = int(x+0.5)
 107
 108 set samples 59
 109 set xrange [1:59]
 110 set auto
 111 set key default
 112
 113 set title "Greatest Common Divisor (for integers only)"
 114
 115 plot gcd(x, 60) with impulses
 116 pause -1 "Hit return to continue"
 117
 118 reset
 119