From b7a2e67ace2cc24ed0beaee1c57e5a095b441d0a Mon Sep 17 00:00:00 2001 From: yaw-man Date: Wed, 8 Feb 2023 15:48:43 -0400 Subject: [PATCH] Feverish insensate spurt: - Added basic AST calculator for derivatives. - Added direct derivation via Faa Di Bruno formula. - Added ansatz method. - Factored out polynomial interpolation functionality. --- ansatz.lua | 21 +++++++++ diff.lua | 110 ++++++++++++++++++++++++++++++++++++++++++++++++ fdb.lua | 103 +++++++++++++++++++++++++++++++++++++++++++++ interpolant.lua | 35 +++++++++++++++ main.lua | 75 ++++++++++++++++++++++++--------- 5 files changed, 325 insertions(+), 19 deletions(-) create mode 100644 ansatz.lua create mode 100644 diff.lua create mode 100644 fdb.lua create mode 100644 interpolant.lua diff --git a/ansatz.lua b/ansatz.lua new file mode 100644 index 0000000..9322695 --- /dev/null +++ b/ansatz.lua @@ -0,0 +1,21 @@ +--Solve the functional equation f' o f = id by ansatz: +--assume f(x) = ax ^ b. Then +--f' o f (x) = ba^b * x ^ ( b ( b-1 ) ) + +local phi = 1 + math.sqrt( 5 ) / 2 +local ihp = 1 - math.sqrt( 5 ) / 2 +local function SeriesCoefficients( ) + + local a, b = math.pow( phi, ihp ), phi + local c, d = a, b + for i = 1, 100 do + + --print( c, d ) + c = c * d + d = d - 1 + + end + +end + +SeriesCoefficients() \ No newline at end of file diff --git a/diff.lua b/diff.lua new file mode 100644 index 0000000..b8c60c9 --- /dev/null +++ b/diff.lua @@ -0,0 +1,110 @@ +--Bespoke symbolic differentiation. +--Little prototypes for building a syntax tree containing only +--multiplication, addition, and composition, +--evaluating its value, and expanding it according to the chain rule. + +local op; +local function d( x ) return x:diff() end +local function D( x ) return setmetatable( { ["n"] = x }, op.atom ) end + +local computedDerivatives = {} + +local __call = D +local __add = function( a, b ) return setmetatable( { a, b }, op.add ) end +local __mul = function( a, b ) return setmetatable( { a, b }, op.mul ) end +local __pow = function( a, b ) return setmetatable( { a, b }, op.compose ) end + +op = { + + atom = { + __add = __add, + __mul = __mul, + __pow = __pow, + __call = __call, + __tostring = function( self ) return tostring( self.n ) end, + + __index = { + + eval = function( self ) + return computedDerivatives[ self.n ] + end, + + diff = function( self ) + return D( self.n + 1 ) + end, + }}, + + compose = { + __add = __add, + __mul = __mul, + __pow = __pow, + __call = __call, + __tostring = function( self ) return "("..tostring( self[1] ).." o "..tostring( self[2] )..")" end, + + __index = { + + eval = function( self ) + --All compositions are of the form f^(n) o f, + --whose value at a fixed point p is f^(n)(p) + return self[1]:eval() + end, + + diff = function( self ) + --All compositions are of the form f^(n) o f, + --whose derivatives are always ( f^(n+1) o f ) + return d( self[1] ) ^ D(0) * D(1) + end, + }}, + + add = { + __add = __add, + __mul = __mul, + __pow = __pow, + __call = __call, + __tostring = function( self ) return tostring( self[1] ).." + "..tostring( self[2] ) end, + + __index = { + + name = "add", + + eval = function( self ) + return self[1]:eval() + self[2]:eval() + end, + + diff = function( self ) + return d(self[1]) + d(self[2]) + end, + }}, + + mul = { + __add = __add, + __mul = __mul, + __pow = __pow, + __call = __call, + __tostring = function( self ) return tostring( self[1] ).." x "..tostring( self[2] ) end, + + __index = { + + eval = function( self ) + return self[1]:eval() * self[2]:eval() + end, + + diff = function( self ) + return d(self[1]) * self[2] + self[1] * d(self[2]) + end, + }}, + + +} + +op.New = function( nodeType, a, b ) + return setmetatable( { a, b }, op[nodeType] or error( "Node Type Not Recognized" ) ) +end + +local a = D(1) ^ D(0) +for i = 1, 5 do + print( a ) + a = a:diff() +end + +return op \ No newline at end of file diff --git a/fdb.lua b/fdb.lua new file mode 100644 index 0000000..d15ba3d --- /dev/null +++ b/fdb.lua @@ -0,0 +1,103 @@ +--Compute derivatives at the fixed point using Faa Di Bruno's formula. + +local ORDER = 21 + +local binCoefs = {} + +local function Choose( n, k ) + if k < 1 then return 1 end + if n < 1 then return 0 end + if k > n then return 0 end + local idx = n * ORDER + k + if binCoefs[idx] then return binCoefs[idx] end + + local z = Choose( n - 1, k - 1) + Choose( n - 1, k ) + binCoefs[ idx ] = z + return z +end + +for i = 13, 1, -1 do + for j = 13, 1, -1 do Choose(i, j) end +end + +Choose = function( n, k ) + local c = binCoefs[ n * ORDER + k ] or ( ( k < 1 ) and 1 or 0 ) + return c +end + +local function Series( p ) + local d = { p, 1.0 / p } -- q[i] := f^(i) (p) + local bellPolynomial = {} + + --Generate incomplete Bell polynomials by recursive formula: + --d[n+1] = - \sum_{k=1}^{n-1} d[k+1] * B(n,k) / p^n + --B[n][k] = \sum_{i=1}^{n-k+1} Choose( n-1, i-1 ) * d[i] * B( n-i, k-1 ) + --sum_i=1^(n-k+1) + + local function Bell( n, k ) + --Base cases. + if n == 0 then + return ( k == 0 ) and 1 or 0 + elseif k == 0 then + return 0 + end + + --Cached values. + local idx = n * ORDER + k + if bellPolynomial[idx] then return bellPolynomial[idx] end + + --Sum. + local b = 0 + for i = 1, n - k + 1 do + local a = Choose( n - 1, i - 1 ) + a = a * d[i] + a = a * Bell( n - i, k - 1 ) + b = b + a + end + bellPolynomial[idx] = b + return b + end + + local NextDerivative = function() + local n = #d + local new = 0 + for k = 1, n - 1 do + new = new + d[k + 1] * Bell(n, k) + end + new = -new / math.pow( p, n ) + d[n + 1] = new + return new + end + + for i = 1, ORDER - 2 do NextDerivative() end + + --Function that evaluates Taylor series. + return { + Interpolant = function( x ) + x = x - p + local pow = 1 + local fact = 1 + local y = p + for i = 1, ORDER do + pow = pow * x + fact = fact * i + y = y + pow * d[i] / fact + end + return y + end, + + Derivative = function( x ) + x = x - p + local pow = 1 + local fact = 1 + local y = 0 + for i = 1, ORDER do + pow = pow * x + fact = fact * i + y = y + pow * d[i] / fact + end + end + } +end + +return Series \ No newline at end of file diff --git a/interpolant.lua b/interpolant.lua new file mode 100644 index 0000000..548d23f --- /dev/null +++ b/interpolant.lua @@ -0,0 +1,35 @@ +--Returns a function which takes a sequence of polynomial coefficients and +--returns two functions which evaluate the polynomial and its derivative, respectively +return function( coefs ) + local fixedPoint = coefs[0] or error( "Must have constant coefficient!" ) + return function(x) + + x = x - fixedPoint + local y = fixedPoint + local pow = 1 + local fact = 1 + for i = 1, #coefs do + pow = pow * x + fact = fact / i + y = y + pow * fact * coefs[i] + end + return y + end, + + function(x) + + x = x - fixedPoint + local y = 0 + local pow = 1 + local fact = 1 + for i = 1, #coefs do + y = y + pow * fact * coefs[i] + pow = pow * x + fact = fact / i + end + return y + end + + + +end diff --git a/main.lua b/main.lua index 6c16bd9..a020bea 100644 --- a/main.lua +++ b/main.lua @@ -7,41 +7,78 @@ --get truncated taylor series expansion of f at p --plot to get some idea about values, convergence. +local RESOLUTION = 1000 +local EXTENT = 16 +local FaaDiBruno = require( 'fdb' ) +local Ansatz = require( 'ansatz' ) +local Poly = require( 'interpolant' ) + +local a = 0.2 local function PlotFunction( f ) - local RESOLUTION = 1000 + local points = {} for n = 1, RESOLUTION do - points[ 2 * n - 1 ] = n / RESOLUTION - points[ 2 * n ] = f( n / RESOLUTION ) + points[ 2 * n - 1 ] = EXTENT * n / RESOLUTION + points[ 2 * n ] = f( EXTENT * n / RESOLUTION ) + + --inverse + --points[ 2 * RESOLUTION + 2 * n ] = points[ 2 * n - 1] + --points[ 2 * RESOLUTION + 2 * n - 1] = points[ 2 * n ] end + local x, y = love.graphics.getDimensions() local tf = love.math.newTransform( - 0, love.graphics.getHeight(), - 0, love.graphics.getWidth(), - -love.graphics.getHeight()) - - love.graphics.setColor( 1, 1, 1, 0.5 ) - love.graphics.setLineWidth( 0.003 ) + x / EXTENT, y,-- / EXTENT, + 0, x / EXTENT, + -y / EXTENT) + local dtf = love.math.newTransform( 0, 0, 0, 1, 1 ) + + --love.graphics.setColor( 1, 1, 1, 0.5 ) + love.graphics.setLineWidth( 0.1 ) love.graphics.setLineJoin( "miter" ) love.graphics.setLineStyle( "smooth" ) local draw = love.draw or function() end love.draw = function() - draw() - love.graphics.replaceTransform( tf ) - love.graphics.line( points ) - end + draw() + love.graphics.push("transform") + love.graphics.translate( 0, y ) + love.graphics.scale( x, -y ) + love.graphics.scale( 1 / EXTENT, 1 / EXTENT ) - love.update = function( dt ) - print( love.mouse.getPosition( ) ) + --love.graphics.scale( 1.0 / EXTENT, 1.0 / EXTENT ) + love.graphics.line( points ) + love.graphics.pop() + + love.graphics.print( a ) end end -PlotFunction( function( x ) return x * x end ) -PlotFunction( function( x ) return x * x * x end ) -PlotFunction( math.sin ) -PlotFunction( function( x ) return math.exp( x ) - 1 end ) \ No newline at end of file +love.update = function( dt ) + if love.keyboard.isScancodeDown("w") then + a = a + 0.05 + love.graphics.setColor( a, 0, 0, a ) + love.draw = nil + --PlotFunction( FaaDiBruno( a ) ) + --PlotFunction( function(x)return x end) + end + + if love.keyboard.isScancodeDown("s") then + a = a - 0.05 + love.graphics.setColor( 1, 0, 0, 1 ) + love.draw = nil + --[[local f = FaaDiBruno(a) + PlotFunction( f.Interpolant ) + PlotFunction( f.Derivative )]] + PlotFunction( function(x)return x end) + local f, df = Poly({[0] = 1.5, 1 + a, 2 }) + PlotFunction( f ) + PlotFunction( df ) + end + if love.keyboard.isScancodeDown("q") then EXTENT = EXTENT + 0.5 end + if love.keyboard.isScancodeDown("e") then EXTENT = EXTENT * 0.99 end +end \ No newline at end of file