diff --git a/ansatz.lua b/ansatz.lua index 9322695..3049102 100644 --- a/ansatz.lua +++ b/ansatz.lua @@ -2,20 +2,31 @@ --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 phi = 0.5 + math.sqrt( 5 ) / 2 +local ihp = 0.5 - math.sqrt( 5 ) / 2 +local a, b = math.pow( phi, ihp ), phi local function SeriesCoefficients( ) local a, b = math.pow( phi, ihp ), phi + local coefs = {[0] = phi} local c, d = a, b - for i = 1, 100 do + for i = 1, 150 do - --print( c, d ) c = c * d d = d - 1 - + coefs[i] = c * math.pow( phi, d ) + end - + return coefs end -SeriesCoefficients() \ No newline at end of file +--SeriesCoefficients() + +local c, d = math.pow( ihp, phi ), ihp +return { + pos = function( x ) return a * math.pow( x, b ) end, + dpo = function( x ) return a * b * math.pow( x, b - 1 ) end, + neg = function( x ) return c * math.pow( x, d ) end, + dne = function( x ) return c * d * math.pow( x, d - 1 ) end, + coefs = SeriesCoefficients() +} \ No newline at end of file diff --git a/fdb.lua b/fdb.lua index d15ba3d..17f1910 100644 --- a/fdb.lua +++ b/fdb.lua @@ -1,6 +1,6 @@ --Compute derivatives at the fixed point using Faa Di Bruno's formula. -local ORDER = 21 +local ORDER = 100 local binCoefs = {} @@ -16,8 +16,8 @@ local function Choose( n, k ) return z end -for i = 13, 1, -1 do - for j = 13, 1, -1 do Choose(i, j) end +for i = ORDER, 1, -1 do + for j = ORDER, 1, -1 do Choose(i, j) end end Choose = function( n, k ) @@ -25,8 +25,10 @@ Choose = function( n, k ) return c end -local function Series( p ) - local d = { p, 1.0 / p } -- q[i] := f^(i) (p) +local function Series( p, order ) + order = order or ORDER + order = math.min( order, ORDER ) + local d = { [0] = p, p, 1.0 / p } -- d[i] := f^(i) (p) local bellPolynomial = {} --Generate incomplete Bell polynomials by recursive formula: @@ -43,7 +45,7 @@ local function Series( p ) end --Cached values. - local idx = n * ORDER + k + local idx = n * order + k if bellPolynomial[idx] then return bellPolynomial[idx] end --Sum. @@ -69,35 +71,11 @@ local function Series( p ) return new end - for i = 1, ORDER - 2 do NextDerivative() 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 - } + --First few derivatives evaluated at fixed point p. + + return d end return Series \ No newline at end of file diff --git a/interpolant.lua b/interpolant.lua index 548d23f..a4c9b62 100644 --- a/interpolant.lua +++ b/interpolant.lua @@ -2,6 +2,7 @@ --returns two functions which evaluate the polynomial and its derivative, respectively return function( coefs ) local fixedPoint = coefs[0] or error( "Must have constant coefficient!" ) + --Interpolant, naive return function(x) x = x - fixedPoint @@ -16,6 +17,7 @@ return function( coefs ) return y end, + --Derivative, naive function(x) x = x - fixedPoint @@ -31,5 +33,4 @@ return function( coefs ) end - end diff --git a/main.lua b/main.lua index a020bea..25c8b9f 100644 --- a/main.lua +++ b/main.lua @@ -7,36 +7,57 @@ --get truncated taylor series expansion of f at p --plot to get some idea about values, convergence. -local RESOLUTION = 1000 -local EXTENT = 16 +local plot = { + x = 0, + X = 1, + y = 0, + Y = 1, + dx = 0.001, + dy = 0.001, + inverse = false, + fdbOrder = 50,} + +plot.Zoom = function( zoomFactor ) + + plot.dx, plot.dy = plot.dx * zoomFactor, plot.dy * zoomFactor + + local cx, cy = 0.5 * ( plot.x + plot.X ), 0.5 * ( plot.y + plot.Y ) + local lx, ly = zoomFactor * 0.5 * ( plot.X - plot.x ), zoomFactor * 0.5 * ( plot.Y - plot.y ) + plot.x, plot.X = cx - lx, cx + lx + plot.y, plot.Y = cy - ly, cy + ly +end + +plot.Translate = function( x, y ) + x = 0.01 * (plot.X - plot.x) * x + y = 0.01 * (plot.Y - plot.y) * y + plot.x, plot.X = plot.x + x, plot.X + x + plot.y, plot.Y = plot.y + y, plot.Y + y +end local FaaDiBruno = require( 'fdb' ) local Ansatz = require( 'ansatz' ) local Poly = require( 'interpolant' ) -local a = 0.2 -local function PlotFunction( f ) - - +local a = 1.619 +local function PlotFunction( f, color ) local points = {} - for n = 1, RESOLUTION do - points[ 2 * n - 1 ] = EXTENT * n / RESOLUTION - points[ 2 * n ] = f( EXTENT * n / RESOLUTION ) + local n = 1 + for x = plot.x - 0.1, plot.X + 0.1, plot.dx do + points[ 2 * n - 1 ] = x + points[ 2 * n ] = f( x ) + n = n + 1 - --inverse - --points[ 2 * RESOLUTION + 2 * n ] = points[ 2 * n - 1] - --points[ 2 * RESOLUTION + 2 * n - 1] = points[ 2 * n ] + end + + --inverse + local inverse = {} + for i = 1, n, 2 do + inverse[i], inverse[i + 1] = points[ i + 1], points[i] end - local x, y = love.graphics.getDimensions() - local tf = love.math.newTransform( - 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 ) + local x, y = love.graphics.getDimensions() love.graphics.setLineWidth( 0.1 ) love.graphics.setLineJoin( "miter" ) love.graphics.setLineStyle( "smooth" ) @@ -47,38 +68,74 @@ local function PlotFunction( f ) love.graphics.push("transform") love.graphics.translate( 0, y ) love.graphics.scale( x, -y ) - love.graphics.scale( 1 / EXTENT, 1 / EXTENT ) + love.graphics.scale( 1 / (plot.X - plot.x), 1/ (plot.Y - plot.y) ) + love.graphics.translate( - plot.x, - plot.y ) + love.graphics.setLineWidth( 0.003 * (plot.X - plot.x) ) - --love.graphics.scale( 1.0 / EXTENT, 1.0 / EXTENT ) - love.graphics.line( points ) + if color then love.graphics.setColor( color ) end + love.graphics.line( plot.inverse and inverse or points ) + love.graphics.circle( "fill", a, a, 0.003 * (plot.X - plot.x) ) love.graphics.pop() love.graphics.print( a ) + love.graphics.print( plot.fdbOrder, 0, 15 ) + love.graphics.print( f(a), 0, 30 ) + love.graphics.print( f(-1.0), 0, 45 ) + love.graphics.setColor( 1,1,1,1 ) + love.graphics.setFont( love.graphics.getFont( 48 ) ) end end -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 +love.wheelmoved = function( x, y ) + return plot.Zoom( (y > 0) and 0.95 or 1.05 ) or love.keypressed() +end + +love.mousepressed = function(x, y, button) + --plot.inverse = not( plot.inverse ) + plot.fdbOrder = plot.fdbOrder + (( button == 1 ) and 1 or -1 ) + return love.keypressed() +end + +--[[local f, df = Poly( FaaDiBruno( a ) ) +PlotFunction( f, {1, 0, 0, 0.5} ) --Function +PlotFunction( df, {0, 1, 0, 0.5} ) --First derivative +PlotFunction( function(x) return x end, {0, 0, 1, 0.5})]] + +love.keypressed = function( key, code ) + if code == "q" then a = a + 0.001 * (plot.X - plot.x) end + if code == "e" then a = a - 0.001 * (plot.X - plot.x) end + if code == "w" then plot.Translate( 0, 1 ) end + if code == "a" then plot.Translate( -1, 0 ) end + if code == "s" then plot.Translate( 0, -1 ) end + if code == "d" then plot.Translate( 1, 0 ) end + if code == "z" then plot.Zoom( 0.95 ) end + if code == "c" then plot.Zoom( 1.05 ) end + if code == "f" then plot.fdbOrder = plot.fdbOrder - 1 end + if code == "r" then plot.fdbOrder = plot.fdbOrder + 1 end + love.draw = nil + + local f, df, fn = Poly( FaaDiBruno( a, plot.fdbOrder ) ) + --PlotFunction( f, {1, 0, 0, 0.3} ) --Function in red. + --PlotFunction( df, {0, 1, 0, 0.3} ) --First derivative in green. + --local lowF, lowDF = Poly( FaaDiBruno( a, 15 ) ) + --PlotFunction( lowF, {1, 0, 0, 0.3} ) --Function in red. + --PlotFunction( lowDF, {0, 1, 0, 0.3} ) --First derivative in green. - 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 ) + --PlotFunction( function(x) return df( f ( x ) ) end, {0, 0, 1, 0.3}) + --PlotFunction( function(x) return f( df ( x ) ) end, {1, 1, 1, 0.3}) + --PlotFunction( function(x) return x end, {1, 1, 1, 0.3}) + local af, an = Poly( Ansatz.coefs ) + PlotFunction( af, {1,1,1, 0.3} ) + PlotFunction( an, {1,1,1, 0.3} ) + PlotFunction( Ansatz.dpo, {1, 1, 1, 0.3} ) + PlotFunction( Ansatz.pos, {1,1,1,0.3} ) +end + +love.update = function( dt ) + for _, code in ipairs{ 'q', 'e', 'w', 'a', 's', 'd', 'z', 'c', 'f', 'r' } do + if love.keyboard.isScancodeDown( code ) then love.keypressed( nil, code ) end 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 +end + +love.keypressed( nil, nil ) \ No newline at end of file