|
Packit |
0986c0 |
#
|
|
Packit |
0986c0 |
# $Id: bivariat.dem,v 1.8 2006/07/06 21:52:19 sfeam Exp $
|
|
Packit |
0986c0 |
#
|
|
Packit |
0986c0 |
# This demo is very slow and requires unusually large stack size.
|
|
Packit |
0986c0 |
# Do not attempt to run this demo under MSDOS.
|
|
Packit |
0986c0 |
#
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
# the function integral_f(x) approximates the integral of f(x) from 0 to x.
|
|
Packit |
0986c0 |
# integral2_f(x,y) approximates the integral from x to y.
|
|
Packit |
0986c0 |
# define f(x) to be any single variable function
|
|
Packit |
0986c0 |
#
|
|
Packit |
0986c0 |
# the integral is calculated using Simpson's rule as
|
|
Packit |
0986c0 |
# ( f(x-delta) + 4*f(x-delta/2) + f(x) )*delta/6
|
|
Packit |
0986c0 |
# repeated x/delta times (from x down to 0)
|
|
Packit |
0986c0 |
#
|
|
Packit |
0986c0 |
delta = 0.2
|
|
Packit |
0986c0 |
# delta can be set to 0.025 for non-MSDOS machines
|
|
Packit |
0986c0 |
#
|
|
Packit |
0986c0 |
# integral_f(x) takes one variable, the upper limit. 0 is the lower limit.
|
|
Packit |
0986c0 |
# calculate the integral of function f(t) from 0 to x
|
|
Packit |
0986c0 |
# choose a step size no larger than delta such that an integral number of
|
|
Packit |
0986c0 |
# steps will cover the range of integration.
|
|
Packit |
0986c0 |
integral_f(x) = (x>0)?int1a(x,x/ceil(x/delta)):-int1b(x,-x/ceil(-x/delta))
|
|
Packit |
0986c0 |
int1a(x,d) = (x<=d*.1) ? 0 : (int1a(x-d,d)+(f(x-d)+4*f(x-d*.5)+f(x))*d/6.)
|
|
Packit |
0986c0 |
int1b(x,d) = (x>=-d*.1) ? 0 : (int1b(x+d,d)+(f(x+d)+4*f(x+d*.5)+f(x))*d/6.)
|
|
Packit |
0986c0 |
#
|
|
Packit |
0986c0 |
# integral2_f(x,y) takes two variables; x is the lower limit, and y the upper.
|
|
Packit |
0986c0 |
# calculate the integral of function f(t) from x to y
|
|
Packit |
0986c0 |
integral2_f(x,y) = (x
|
|
Packit |
0986c0 |
-int2(y,x,(x-y)/ceil((x-y)/delta))
|
|
Packit |
0986c0 |
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.)
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
set autoscale
|
|
Packit |
0986c0 |
set title "approximate the integral of functions"
|
|
Packit |
0986c0 |
set samples 50
|
|
Packit |
0986c0 |
set key bottom right
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
f(x) = exp(-x**2)
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
plot [-5:5] f(x) title "f(x)=exp(-x**2)", \
|
|
Packit |
0986c0 |
2/sqrt(pi)*integral_f(x) title "erf(x)=2/sqrt(pi)*integral_f(x)", \
|
|
Packit |
0986c0 |
erf(x) with points
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
pause -1 "Hit return to continue"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
f(x)=cos(x)
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
plot [-5:5] f(x) title "f(x)=cos(x)", integral_f(x)
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
pause -1 "Hit return to continue"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
set title "approximate the integral of functions (upper and lower limits)"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
f(x)=(x-2)**2-20
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
plot [-10:10] f(x) title "f(x)=(x-2)**2-20", integral2_f(-5,x)
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
pause -1 "Hit return to continue"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
f(x)=sin(x-1)-.75*sin(2*x-1)+(x**2)/8-5
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
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)
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
pause -1 "Hit return to continue"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
#
|
|
Packit |
0986c0 |
# This definition computes the ackermann. Do not attempt to compute its
|
|
Packit |
0986c0 |
# values for non integral values. In addition, do not attempt to compute
|
|
Packit |
0986c0 |
# its beyond m = 3, unless you want to wait really long time.
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
ack(m,n) = (m == 0) ? n + 1 : (n == 0) ? ack(m-1,1) : ack(m-1,ack(m,n-1))
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
set xrange [0:3]
|
|
Packit |
0986c0 |
set yrange [0:3]
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
set isosamples 4
|
|
Packit |
0986c0 |
set samples 4
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
set title "Plot of the ackermann function"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
splot ack(x, y)
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
pause -1 "Hit return to continue"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
set xrange [-5:5]
|
|
Packit |
0986c0 |
set yrange [-10:10]
|
|
Packit |
0986c0 |
set isosamples 10
|
|
Packit |
0986c0 |
set samples 100
|
|
Packit |
0986c0 |
set key top right at 4,-3
|
|
Packit |
0986c0 |
set title "Min(x,y) and Max(x,y)"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
#
|
|
Packit |
0986c0 |
min(x,y) = (x < y) ? x : y
|
|
Packit |
0986c0 |
max(x,y) = (x > y) ? x : y
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
plot sin(x), x**2, x**3, max(sin(x), min(x**2, x**3))+0.5
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
pause -1 "Hit return to continue"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
#
|
|
Packit |
0986c0 |
# gcd(x,y) finds the greatest common divisor of x and y,
|
|
Packit |
0986c0 |
# using Euclid's algorithm
|
|
Packit |
0986c0 |
# as this is defined only for integers, first round to the nearest integer
|
|
Packit |
0986c0 |
gcd(x,y) = gcd1(rnd(max(x,y)),rnd(min(x,y)))
|
|
Packit |
0986c0 |
gcd1(x,y) = (y == 0) ? x : gcd1(y, x - x/y * y)
|
|
Packit |
0986c0 |
rnd(x) = int(x+0.5)
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
set samples 59
|
|
Packit |
0986c0 |
set xrange [1:59]
|
|
Packit |
0986c0 |
set auto
|
|
Packit |
0986c0 |
set key default
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
set title "Greatest Common Divisor (for integers only)"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
plot gcd(x, 60) with impulses
|
|
Packit |
0986c0 |
pause -1 "Hit return to continue"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
#
|
|
Packit |
0986c0 |
# This definition computes the sum of the first 10, 100, 1000 fourier
|
|
Packit |
0986c0 |
# coefficients of a (particular) square wave.
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
set title "Finite summation of 10, 100, 1000 fourier coefficients"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
set samples 500
|
|
Packit |
0986c0 |
set xrange [-10:10]
|
|
Packit |
0986c0 |
set yrange [-0.4:1.2]
|
|
Packit |
0986c0 |
set key bottom right
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
fourier(k, x) = sin(3./2*k)/k * 2./3*cos(k*x)
|
|
Packit |
0986c0 |
sum10(x) = 1./2 + sum [k=1:10] fourier(k, x)
|
|
Packit |
0986c0 |
sum100(x) = 1./2 + sum [k=1:100] fourier(k, x)
|
|
Packit |
0986c0 |
sum1000(x) = 1./2 + sum [k=1:1000] fourier(k, x)
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
plot \
|
|
Packit |
0986c0 |
sum10(x) title "1./2 + sum [k=1:10] sin(3./2*k)/k * 2./3*cos(k*x)", \
|
|
Packit |
0986c0 |
sum100(x) title "1./2 + sum [k=1:100] sin(3./2*k)/k * 2./3*cos(k*x)", \
|
|
Packit |
0986c0 |
sum1000(x) title "1./2 + sum [k=1:1000] sin(3./2*k)/k * 2./3*cos(k*x)"
|
|
Packit |
0986c0 |
pause -1 "Hit return to continue"
|
|
Packit |
0986c0 |
|
|
Packit |
0986c0 |
reset
|
|
Packit |
0986c0 |
|