Large commit, fixed large mistake in calculating Ansatz values, added tests for FDB series plot.
This commit is contained in:
parent
b7a2e67ace
commit
98c27e0725
23
ansatz.lua
23
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()
|
||||
--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()
|
||||
}
|
||||
44
fdb.lua
44
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,
|
||||
--First few derivatives evaluated at fixed point p.
|
||||
|
||||
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
|
||||
}
|
||||
return d
|
||||
end
|
||||
|
||||
return Series
|
||||
|
|
@ -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
|
||||
|
|
|
|||
143
main.lua
143
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
|
||||
|
||||
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 )
|
||||
--inverse
|
||||
local inverse = {}
|
||||
for i = 1, n, 2 do
|
||||
inverse[i], inverse[i + 1] = points[ i + 1], points[i]
|
||||
end
|
||||
|
||||
--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.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.
|
||||
|
||||
--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 )
|
||||
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)
|
||||
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
|
||||
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
|
||||
love.keypressed( nil, nil )
|
||||
Loading…
Reference in New Issue