1
0
Fork 0
Libraries for bc and dc.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

112 lines
3.0 KiB

#!/usr/local/bin/bc -l
### IntDiff.BC - numeric differentiation and integration of
### a single variable function
define f(x) { return x^2 } # example - redefine in later code/bc session
# Numerically differentiate the global function f(x)
define dfxdx(x) {
auto eps;
eps = A^-scale
scale *= 2
x = (f(x+eps)-f(x-eps))/(2*eps)
scale /= 2
return x/1
}
# New global variable like 'scale' - determines accuracy of numerical
# integration at the expense of time. Don't set above 15 unless this
# is running on a really fast machine!
depth = 10
# Numerically integrate the global function f(x) between x = a and x = b
# . uses the trapezoidal rule
define ifxdx_old(a,b) {
auto os,h,s,t,i
if(a==b)return f(a)
os = scale;if(scale<depth)scale=depth
scale+=3
h = 2^-depth
if(b<a){i=b;b=a;a=i}
s = (b - a) * h
t =(f(a)+f(b))/2
for(i=a+s;i<b;i+=s)t+=f(i)
scale=os;return t*s/1
}
# Numerically integrate the global function f(x) between x = a and x = b
define ifxdx(a,b) {
auto oib,od,os,s,s8,t,i,j,ni,fi,fis
if(a==b)return f(a)
od=depth;if(depth<3)depth=3
os=scale;if(scale<(i=depth+depth))scale=i
scale+=3
if(b<a){i=b;b=a;a=i}
s=(b-a)*(2^-depth)
oib=ibase;ibase=A
s8 = s*8
fi = 989*f(a)
for(i=a;i<b;i=ni){
ni=j=i+s8;
t+= fi+(fis=989*f(j))
t+= 5888*(f(i+=s)+f(j-=s))
t-= 928*(f(i+=s)+f(j-=s))
t+=10496*(f(i+=s)+f(j-=s))
t-= 4540*f(i+=s)
fi=fis
}
depth=od;scale=os
t*=s*4/14175
ibase=oib;return t
}
# glai - guess limit at infinity
# Assumes p, q and r are 3 consecutive convergents to a limit and
# attempts to extrapolate precisely what that limit is after an infinite
# number of iterations.
# 0 = glai returns function result only
# 1 = glai commentates on interesting convergents
glaitalk = 1
define glai(p,q,r) {
auto m,n
m = q^2-p*r
n = 2*q-p-r
if(n==0)if(m==0){
if(glaitalk)print "glai: Constant series detected\n"
return p
}else{
if(glaitalk)print "glai: Arithmetic progression detected: limit is infinite\n"
return 1/0
}
if(m==0){
if(glaitalk)print "glai: Geometric progression detected: limit wraps to zero!\n"
return 0
}
return m/n
}
# Examples:
# glai(x,x+1,x+2) causes a division by zero error as the limit of
# an arithmetic progression is infinite
# glai(a*k,a*k^2,a*k^3) returns zero! The limit of a geometric
# progression in p-adics is precisely that,
# and somehow this simple function 'knows'.
# glai(63.9, 63.99, 63.999) returns 64 - correctly predicting the
# limit of the sequence.
# Run consecutive convergents to the ifxdx function through glai()
# attaining "better" accuracy with slightly fewer calculations
define ifxdx_g(a,b) {
auto p,q,r
depth-=3 ; p = ifxdx(a,b)
.=depth++ ; q = ifxdx(a,b)
.=depth++ ; r = ifxdx(a,b)
.=depth++
return glai(p,q,r)
}
#define f(x){if(x<=0)return 0;x=root(x,x);return x*(x-1)}
#zz=-0.10717762842559665710112408473270837028206726160094438