juliaapproximation / approxfun.jl Goto Github PK
View Code? Open in Web Editor NEWJulia package for function approximation
Home Page: http://juliaapproximation.github.io/ApproxFun.jl/
License: Other
Julia package for function approximation
Home Page: http://juliaapproximation.github.io/ApproxFun.jl/
License: Other
Is something like this possible?
d = Interval(0.,50.)
D = diff(d)
t = Fun(identity,d)
F = D^2 +.5D + 1
BC = [ldirichlet(d), lneumann(d)]
BC0 = [1.0, 0.0]
x = [BC, F]\[BC0, 0.0] #evolution of undriven equation x'' + .5x' + x = 0, works beautifully
A= [ ldirichlet(d) 0;
lneumann(d) 0;
rdirichlet(d) 0;
F -1; ]
b = [1.0; 0.0; 2.0; 0.0]
xu = A\b
xu
should be two Funs corresponding to x'' + .5x' + x = u
with the desired boundary conditions.
Does that seem reasonable?
Is there a reason why these shouldn't be identical operators
julia> evaluate(d,d.a)
Evaluation{UltrasphericalSpace{0},Float64,Float64}(UltrasphericalSpace{0}(Interval{Float64}(0.0,50.0)),0.0,0)
julia> ldirichlet(d)
Evaluation{UltrasphericalSpace{0},Bool,Float64}(UltrasphericalSpace{0}(Interval{Float64}(0.0,50.0)),false,0)
for any d, in this case d=Interval(0.,50.)
At least they should behave the same, but something is up here
d = Interval(0.,50.)
D = diff(d)
t = Fun(identity,d)
F = D^2 +.5D + 1
BC_test1 = [evaluate(d,0.), evaluate(d,0.)*D]
BC_test2 = [ldirichlet(d), lneumann(d)]
BC0 = [1.0, 0.0]
x1 = [BC_test1, F]\[BC0, 0.0]
x2 = [BC_test2, F]\[BC0, 0.0]
julia> norm(x1-x2)
0.9875181589473826
ldirichlet(d)
and evaluate(d,d.a)
behave differently and/or lneumann(d)
and evaluate(d,d.a) * D
I know the PiecewiseSpace isn't exported yet, but I was experimenting with it a little bit. Consider a simple absolute value override:
function Base.abs(f::Fun{ChebyshevSpace})
d=domain(f)
if length(f) <=1
pts=[d.a,d.b]
else
pts=roots(f)
isapprox(d.a,pts[1]) ? pts[1] = d.a : pts = [d.a,pts]
isapprox(d.b,pts[end]) ? pts[end] = d.b : pts = [pts,d.b]
end
Fun(x->abs(f[x]),pts)
end
Then the following code doesn't work for a couple different reasons:
x = Fun(identity)
sf(i) = sinpi(i.x)
absf(i) = abs(sf(i))
for i=1:9
println(absf(i)[.51]," ",abs(sinpi(i..51)))
end
There seems to be two issues: the first occurring at i=5 in the evaluation, perhaps due to erroneous interlacing of the coefficients, and the second at the end, where it crashes. A third issue is passing the Fun an array argument:
xx = linspace(-1,1,11);
absf(1)[xx]
julia> s = Fun(sin, [0., 7.]);
julia> c = Fun(cos, [0., 7.]);
julia> as = abs(s);
julia> ac = abs(c);
julia> ac+as;
ERROR: assertion failed: domainscompatible(a,b)
in conversion_rule at /Users/goretkin/.julia/v0.3/ApproxFun/src/Spaces/Ultraspherical/UltrasphericalOperators.jl:273
in maxspace at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/FunctionSpace.jl:109
in maxspace at /Users/goretkin/.julia/v0.3/ApproxFun/src/Spaces/Modifier/VectorSpace.jl:149
in + at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/Fun.jl:119
I spent a couple of hours trying to figure out how I should implement this but I don't know enough about ApproxFun. I put what I started in a gist https://gist.github.com/goretkin/5d4fb776b9ef09f4bd1f
I think the point here is to add ac + as
without recomputing the coefficients, though one stopgap could be to override the addition method on Fun{PiecewiseSpace}
directly and compute Fun( (x)-> ac[x] + as[x], unique(merge_sorted(breaks(ac.space), breaks(as.space))))
Probably due to your recent refactoring, I'm unable to plot a Fun
now:
f = Fun(x->exp(x.^2),[-10,10])
plot(f)
`plot` has no method matching plot(::Fun{UltrasphericalSpace{0},Float64})
while loading In[11], in expression starting on line 2
I did find some anomalous behaviour, though. Consider for instance:
d = Interval()^2
x = Fun((x,y)->x,d)
y = Fun((x,y)->y,d)
The x TensorFun is fine, but y fails to instantiate.
Maybe there's already a way to do this and I don't see.
Say x
is some Fun. I would like to project it into the space of functions f | f[0] = 0
. That's a linear subspace. There needs to be some norm to define "project".
I was about to construct a projection matrix for the coefficients of a given length, but maybe I can do this in terms of existing machinery.
This would allow a Dirichlet functional, which adds continuity conditions for piecewise spaces
Does this package need to be dependent on PyPlot or any plotting package?
And another one is:
T(n) = Fun(x->cos(n.*acos(x)))
T(7)
T(8)
I believe this has to do with the initial number of sample points taken in the constructor and the Chebyshev aliasing. I think the Chebfun solutions is to add one rand() sample check on the interval to verify if it is ok to stop looping?
The convention for Julia packages is that the repositories end in .jl
.
d = Interval(0.,50.)
D = diff(d)
t = Fun(identity,d)
F = D^2 +.5D + 1
BC = [ldirichlet(d), lneumann(d)]
BC0 = [1.0, 0.0]
x = [BC, F]\[BC0, (sin(3t))] #evolution of driven equation x'' + .5x' + x = sin(3t)
doesn't work if the driving term is piecewise.
julia> x = [BC, F][BC0, abs(sin(3t))]
ERROR: points
has no method matching points(::UnionDomain{Interval{Float64}}, ::Int64, ::Type{Float64})
in Fun at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/constructors.jl:4
in zerocfsFun at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/constructors.jl:97
in Fun at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/constructors.jl:143
in Fun at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/constructors.jl:154
in Fun at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/constructors.jl:47
in Fun at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/constructors.jl:177
in abs at /Users/goretkin/.julia/v0.3/ApproxFun/src/Extras/specialfunctions.jl:32
Probably needs to wait til 0.4
Fun(x->x<0?0:1,[-1.,0.,1.])
breaks cause second kind points evaluate at endpoints. I can't see a reason not to use first kind points which avoid endpoints.
Simple example of under-determined system
d = Interval()
[ldirichlet(d),diff(d)] \ [0, zero] #under-determined
This was the example of the over-determined system that prompted me posting the issue
d = Interval()
D = diff(d)
F = D^2 +.5D + 1
A= [ ldirichlet(d) 0;
lneumann(d) 0;
rdirichlet(d) 0;
ldirichlet(d) 0;
F -1; ]
b = [1.0; 0.0; 2.0; 0.0; 0.0]
julia> xu = A\b
WARNING: Maximum number of iterations 1000000 reached
Fun{VectorDomainSpace{2,UltrasphericalSpace{0},Float64},Float64}([NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN … NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN],VectorDomainSpace{2,UltrasphericalSpace{0},Float64}(UltrasphericalSpace{0}(Interval{Float64}(-1.0,1.0))))
(by the way, I had meant to have)
A= [ ldirichlet(d) 0;
lneumann(d) 0;
rdirichlet(d) 0;
rneumann(d) 0;
F -1; ]
I wish I could type
D=Derivative()
B=dirichlet()
x=Fun(identity)
u=[B,D^2+x][1.,0.]
and it chooses the space from x. At the moment, the x is interpreted as Multiplication(x,AnySpace()) so that spaces can be inferred from other operators. (Since Chebyshev can multiply other spaces like Jacobi). Some system for choosing a default space should exist when it returns AnySpace.
not sure if this is a legitimate issue, but
julia> Fun(x->2,d)
ERROR: type: chebyshevtransform: in typeassert, expected Array{Int64,1}, got Array{Float64,1}
in chebyshevtransform at /Users/goretkin/.julia/v0.3/ApproxFun/src/LinearAlgebra/chebyshevtransform.jl:32
in chebyshevtransform at /Users/goretkin/.julia/v0.3/ApproxFun/src/LinearAlgebra/chebyshevtransform.jl:27
in transform at /Users/goretkin/.julia/v0.3/ApproxFun/src/Spaces/Ultraspherical/ChebyshevSpace.jl:22
in Fun at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/constructors.jl:10
in zerocfsFun at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/constructors.jl:89
in Fun at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/constructors.jl:137
in Fun at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/constructors.jl:142
The error is not directly related to it, but it might trip someone up that they see expected Array{Int64,1}, got Array{Float64,1}
, when in fact the issue is that the function is returning Int64 instead of Float64
The citation
S. Olver & A. Townsend (2013), A fast and well-conditioned spectral method, SIAM Review, 55:462Ð489
in the README.md contains the strange character Ð
. This should probably be a hyphen?
Hi Sheehan,
I'm trying to create the Volterra integral operator and am running into difficulties. Since Integral(d) is an indefinite antiderivative and the Volterra integral operator is from -1...x, it ends up being a banded operator with a boundary row and may need an algebra BandedOperator + Functional....
For instance, suppose we want to solve \int_{-1}^x u(y) dy + u(x) = 1, with solution u(x) = e^{-1-x}.
using ApproxFun
x = Fun(identity)
d = domain(x)
f = Fun([1.0],d)
Volt = Integral(d)
L = Volt + I
u = L\f
does not give the correct solution. Instead,
k = 20
Lp = full(L[1:k,1:k]) - [(ldirichlet(d)*Volt)[1:k]';zeros(k-1,k)]
u = Fun(Lp\[pad(f.coefficients,20)],UltrasphericalSpace{1}(d))
will give the correct solution. I think some kind of algebra is needed to add Functionals and BandedOperators, because now Lp could be represented as an Array{Operator,1} with a Functional and a BandedOperator. However, the addentries! would have to know to shift the row index by 1.
How exactly should someone use plan_chebyshevtransform()?
I get the following:
v = rand(10)
exact = chebyshevtransform(v)
p = plan_chebyshevtransform(v)
x1 = chebyshevtransform(v,p)
x2 = p(v)
println(norm(x1 - exact))
println(norm(x2 - exact))
0.0
11.423085769164677
From our email exchange, perhaps you expected p(v)
to work and there is a bug?
clenshaw
tries to take imag
of non-real numbers, and breaks evaluation at a complex x
into
clenshaw(c,real(x))+im*clenshaw(c,imag(x))
. Does this correspond to evaluating the polynomial at a complex point or is it something else?
I am trying to evaluate the polynomial at a Dual number, and I don't think this is the behavior I want.
When there is no ambiguity, I think this makes things cleaner
Some operators have definitions for any spaces, e.g. ConstantOperator is always diagonal. These have range/domain space as AnySpace. But other operators are not defined until spaces are attached, e.g. Derivative. At the moment they default to returning AnySpace (until a space is inferred by say additition), but really they should be returning something like UnknownSpace.
currently gives a stack overflow
julia> x = Fun(identity)
Fun{UltrasphericalSpace{0},Float64}([0.0,1.0],UltrasphericalSpace{0}(Interval{Float64}(-1.0,1.0)))
julia> x[1.0:0.1:5.0]
ERROR: stack overflow
If you evaluate the FloatRange
as an array before ( [1.:.1:5.]
), then there's no problem.
This stack trace I wasn't able to reproduce may help
julia> vec(yn)[1][3:.01:4]
ERROR: stack overflow
in abs at operators.jl:359
in clenshaw at /Users/goretkin/.julia/v0.3/ApproxFun/src/LinearAlgebra/clenshaw.jl:12
in clenshaw at /Users/goretkin/.julia/v0.3/ApproxFun/src/LinearAlgebra/clenshaw.jl:15 (repeats 52398 times)
in getindex at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/Fun.jl:76
matplotlib plotting is very slow. May GLPlot is better?
for consideration,
diff of a vector-domain function should probably raise a warning or be redefined to be
as approximately
diff{T<:VectorDomainSpace} (x::Array{Fun{T,...}}) = diff(devec(x)
Right now the space promotion system consists of two parallel but intertwined systems, which might be too complicated.
A) Converting coefficients (no need for banded operator)
spaceconversion(cfs,A,B): converts between coefficients in space A to coefficients space B
unionspace(A,B): return a space that both space A and space B can be converted
B) Constructing banded Conversion operator
Conversion(A,B): construct banded conversion operator
maxspace(A,B): return space that has banded conversion operator from A and B
minspace(A,B): return space that has banded conversion operator to A and B
conversion_type(A,B)/conversion_rule(A,B): simplifies implementation of maxspace/minspace when the returned type is either A or B
canonicalspace(A): returns a space that is known there exists banded conversion to/from that can be used to find out how to convert to other space (usually ChebyshevSpace)
julia> d = Interval()
Interval{Float64}(-1.0,1.0)
julia> d = Interval();
julia> t = Fun(identity,d);
julia> f = devec([t^2, sin(t)]);
julia> diff(d) * f #should be the same as diff(f)
ERROR: `canonicalspace` has no method matching canonicalspace(::VectorDomainSpace{2,UltrasphericalSpace{0},Float64})
in spaceconversion at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/FunctionSpace.jl:185
in * at /Users/goretkin/.julia/v0.3/ApproxFun/src/Operators/algebra.jl:400
julia> eye(2) * f #should be the same as eye(2) * vec(f)
ERROR: `*` has no method matching *(::Array{Float64,2}, ::Fun{VectorDomainSpace{2,UltrasphericalSpace{0},Float64},Float64})
and I guess we would have to allow
julia> eye(2) * diff(d) #and diff(d) * eye(2)
ERROR: `*` has no method matching *(::Array{Float64,2}, ::Derivative{UltrasphericalSpace{0},Float64})
and somehow allow the resulting operator to only apply to vector-valued funs of the correct dimension.
I recently installed Julia on my institution supplied desktop, and after adding ApproxFun and loading it, I get this error:
ERROR: T not defined
in include at ./boot.jl:242
in include_from_node1 at ./loading.jl:128
in include at ./boot.jl:242
in include_from_node1 at ./loading.jl:128
in include at ./boot.jl:242
in include_from_node1 at ./loading.jl:128
in reload_path at ./loading.jl:152
in _require at ./loading.jl:67
in require at ./loading.jl:52
in require_3B_3936 at /home/slevinsky/Documents/julia/usr/bin/../lib/julia/sys.so
while loading /home/slevinsky/.julia/v0.4/ApproxFun/src/Operators/spacepromotion.jl, in expression starting on line 3
while loading /home/slevinsky/.julia/v0.4/ApproxFun/src/Operators/Operator.jl, in expression starting on line 116
while loading /home/slevinsky/.julia/v0.4/ApproxFun/src/ApproxFun.jl, in expression starting on line 22
It seems that the problem is related to the type declarations of the operators. In the current build, type parameters cannot take other parameters as parameters. I believe JuliaArrays/ArrayViews.jl#15 is a similar issue, but don't know if the same fix applies or what impact it will have on ApproxFun.
With the current Julia master, I get:
julia> using ApproxFun
ERROR: T not defined
in include at ./boot.jl:242
in include_from_node1 at ./loading.jl:128
in include at ./boot.jl:242
in include_from_node1 at ./loading.jl:128
in include at ./boot.jl:242
in include_from_node1 at ./loading.jl:128
in reload_path at ./loading.jl:152
in _require at ./loading.jl:67
in require at ./loading.jl:52
in require at /Users/viral/julia/usr/lib/julia/sys.dylib
while loading /Users/viral/.julia/v0.4/ApproxFun/src/Operators/spacepromotion.jl, in expression starting on line 3
while loading /Users/viral/.julia/v0.4/ApproxFun/src/Operators/Operator.jl, in expression starting on line 116
while loading /Users/viral/.julia/v0.4/ApproxFun/src/ApproxFun.jl, in expression starting on line 22
I keep getting the same error, either in Terminal or in IJulia. Has anyone else run into this? I'm on a Macbook Pro OS X 10.9, Julia Version 0.4.0-dev+732 Commit 5d8c6d3*, all packages up to date.
using GLPlot
using ApproxFun
d = Interval()^2
u0 = Fun((x,y)->exp(-40(x-.1)^2-40(y+.2)^2),d)
B = dirichlet(d)
D = Derivative(Interval())
L = (0.01D^2-4D)⊗I + I⊗(0.01D^2-3D)
h = 0.002
timeevolution(B,L,u0,h)
ERROR: key not found: Window(GLPlot)
in update at /Users/Mikael/.julia/v0.4/GLWindow/src/reactglfw.jl:69
in window_position at /Users/Mikael/.julia/v0.4/GLWindow/src/reactglfw.jl:90
in createwindow at /Users/Mikael/.julia/v0.4/GLWindow/src/reactglfw.jl:264
in createdisplay at /Users/Mikael/.julia/v0.4/GLPlot/src/GLPlot.jl:46
in glsurf at /Users/Mikael/.julia/v0.4/ApproxFun/src/Plot/GLPlot.jl:57
in surf at /Users/Mikael/.julia/v0.4/ApproxFun/src/Plot/Plot.jl:66
in timeevolution at /Users/Mikael/.julia/v0.4/ApproxFun/src/Extras/timeevolution.jl:137
in timeevolution at /Users/Mikael/.julia/v0.4/ApproxFun/src/Extras/timeevolution.jl:140
Right now a constant is f=IFun([c],d). But points(f) returns NaN. Should it be 0?
The following seems like normal and expected behavior. Maybe it could be used in the documentation to illustrate that truncating the coefficients isn't just a matter of efficiency.
julia> f = Fun(sin,[0.0,7.0])
Fun{UltrasphericalSpace{0},Float64}([0.133342,-0.257296,0.321759,0.724387,-0.143404,-0.150661,0.0178401,0.0126291,-0.00108255,-0.00058234 … -9.52795e-7,-3.48357e-7,1.65558e-8,5.22381e-9,-2.16483e-10,-6.00997e-11,2.20858e-12,5.47562e-13,-1.81799e-14,-3.96905e-15],UltrasphericalSpace{0}(Interval{Float64}(0.0,7.0)))
julia> plot(f)
julia> plot(diff(f))
julia> plot(diff(f,20))
Understandably, after taking many many derivatives of a degree-22 polynomial, the derivative eventually vanishes.
If the function is instead constructed as
f = Fun(sin,[0.0,7.0],200)
, differentiation seems to amplify the noise. plot(diff(f,5))
is very wrong.
The same type of thing happens with the Fourier basis
julia> f = FFun(sin,200)
Fun{PeriodicSumSpace{Complex{Float64},HardySpace{true},HardySpace{false}},Complex{Float64}}(Complex{Float64}[-6.12323e-19+0.0im,-2.0345e-17+0.5im,-1.24328e-17-0.5im,-6.12323e-19-2.03847e-17im,-6.12323e-19+2.03847e-17im,5.80908e-18+7.24421e-18im,5.80908e-18-7.24421e-18im,-6.12323e-19-1.02324e-17im,-6.12323e-19+1.02324e-17im,7.38566e-19+2.0287e-18im … 2.5834e-18+1.00975e-17im,-6.12323e-19+2.61691e-17im,-6.12323e-19-2.61691e-17im,-1.13223e-18-2.89214e-17im,-1.13223e-18+2.89214e-17im,-6.12323e-19+7.03922e-18im,-6.12323e-19-7.03922e-18im,2.55521e-18-0.0im,9.01674e-18-1.98692e-17im,-6.12323e-19+0.0im],PeriodicSumSpace{Complex{Float64},HardySpace{true},HardySpace{false}}((HardySpace{true}(PeriodicInterval{Float64}(-3.141592653589793,3.141592653589793)),HardySpace{false}(PeriodicInterval{Float64}(-3.141592653589793,3.141592653589793)))))
julia> plot(diff(f,7))
I tried to implement the gamma function from its integral definition but it seemed to choke:
Gamma(x) = cumsum(Fun(t->t^(x-1)*exp(-t),[0.0,Inf]))[Inf]
Gamma(5)
WARNING: Maximum length reached
ERROR: `chebyshevintegrate` has no method matching chebyshevintegrate(::Ray, ::Array{Float64,1})
in cumsum at /Users/jiahao/.julia/v0.3/ApproxFun/src/Fun/Fun.jl:185
in Gamma at none:1
Gamma100(x) = cumsum(Fun(t->t^(x-1)*exp(-t),[0.0,100.0]))[100.0] #Approximate true gamma function
@show [Gamma100(5) factorial(4)]
@show [Gamma100(0.5) π/2]
[Gamma100(5) factorial(4)] => [24.000000000000014 24.0]
WARNING: Maximum length reached
[Gamma100(0.5) π / 2] => [NaN 1.5707963267948966]
In the Readme, regarding PDEs:
d=Interval()^2
.
In Julia 0.3.0, ApproxFun 0.0.3, This raises an error,
ERROR: '*' has no method matching *(::Interval{Float64}, ::Interval{Float64}) in power_by_squaring at intfuncs.jl:56 in ^ at intfuncs.jl:86
And somewhat hobbles the PDE solving efforts. Indeed, it appears no arithmetic is defined for the Interval
type. Is there an implicit dependency that one must build?
The current approach for multiple pieces is slow. It can be sped up with the following iterative approach.
function schursolve(Lin,rhs)
@Assert length(rhs)<=2 # TODO: forcing term
if !isa(domainspace(Lin),PiecewiseSpace)
return [dirichlet(domainspace(Lin)),Lin]\rhs
end
L=blkdiag(Lin)
if length(L)==2
L1,L2=L[1],L[2]
elseif length(L)==3
L1,L2=L[1],blkdiagm(L[2:3])
else #length>3
d=div(length(L),2)
L1,L2=blkdiagm(L[1:d]),blkdiagm(L[d+1:end])
end
sp1=domainspace(L1);sp2=domainspace(L2)
c1=rhs[1];
c2=length(rhs)>=2?rhs[2]:0.
u,v=neumann(sp1)[2],-neumann(sp2)[1]
f=[0.,1.];g=[1.] # these are the patching conditions
Af=schursolve(L1,f); Bg=schursolve(L2,g )
uAf=u*Af; vBg=v*Bg
Σ=-uAf-vBg
Ac1=schursolve(L1,[c1]);Bc2=schursolve(L2,[0.,c2])
uAc1=u*Ac1;vBc2=v*Bc2
UΓ=(uAc1+vBc2)/Σ
U1=Ac1+UΓ*Af
U2=Bc2+UΓ*Bg
U=devec([isa(space(U1),PiecewiseSpace)?vec(U1):U1,
isa(space(U2),PiecewiseSpace)?vec(U2):U2])
end
I don't know exactly how to ask this question, but I'd like to try some optimization over Fun
s.
I found this paper: http://www-m3.ma.tum.de/foswiki/pub/M3/ChebOpt/WebHome/Minimization_of_Variational_Integrals_using_the_Chebfun_System.pdf but I haven't read through it.
I basically would like to try out taking some first-order constrained optimization algorithm (maybe a penalty method) and replace the decision variable vector with an ApproxFun. Does that make sense?
julia> x = Fun(identity,Interval(0.,1.))
Fun{UltrasphericalSpace{0},Float64}([0.5,0.5],UltrasphericalSpace{0}(Interval{Float64}(0.0,1.0)))
julia> @which vec(abs(sin(x)))
vec{S<:DomainSpace{T},V,T}(f::Fun{PiecewiseSpace{S<:DomainSpace{T},V},T}) at /Users/goretkin/.julia/v0.3/ApproxFun/src/Spaces/Modifier/VectorSpace.jl:56
julia> @which vec(ApproxFun.devec([sin(x),cos(x)]))
vec{n,S<:DomainSpace{T},V,T}(f::Fun{VectorDomainSpace{n,S<:DomainSpace{T},V},T}) at /Users/goretkin/.julia/v0.3/ApproxFun/src/Spaces/Modifier/VectorSpace.jl:26
For a piecewise Fun, vec
returns the individual pieces as separate Funs. For a vector Fun, it returns an Array of Funs.
Not sure what is best here, but I think the first should be something like pieces
. I don't know if pieces
should work on a vector-valued Fun, where each entry is Piecewise (the boundaries of each entry might not align. These "modifiers" might not interact nicely right now.
julia> vecpiece = ApproxFun.devec([abs(sin(x)),abs(cos(x))])
Fun{PiecewiseSpace{PiecewiseSpace{UltrasphericalSpace{0},Float64},Float64},Float64}([1.33227e-15,0.602195,0.472001,0.472001,0.472001,0.472001,0.472001,0.643363,0.277826,-0.513625 … -2.92839e-16,0.0,0.0,0.0,2.00352e-15,0.0,1.9792e-15,0.0,1.92415e-15,0.0],PiecewiseSpace{PiecewiseSpace{UltrasphericalSpace{0},Float64},Float64}([PiecewiseSpace{UltrasphericalSpace{0},Float64}([UltrasphericalSpace{0}(Interval{Float64}(0.0,2.6645352591003757e-15)),UltrasphericalSpace{0}(Interval{Float64}(2.6645352591003757e-15,3.1415926535897922)),UltrasphericalSpace{0}(Interval{Float64}(3.1415926535897922,6.283185307179586)),UltrasphericalSpace{0}(Interval{Float64}(6.283185307179586,9.42477796076938)),UltrasphericalSpace{0}(Interval{Float64}(9.42477796076938,10.0))]),PiecewiseSpace{UltrasphericalSpace{0},Float64}([UltrasphericalSpace{0}(Interval{Float64}(0.0,1.5707963267948983)),UltrasphericalSpace{0}(Interval{Float64}(1.5707963267948983,4.712388980384691)),UltrasphericalSpace{0}(Interval{Float64}(4.712388980384691,7.853981633974486)),UltrasphericalSpace{0}(Interval{Float64}(7.853981633974486,10.0))])]))
julia> vecpiece[2]
ERROR: `start` has no method matching start(::UnionDomain{Interval{Float64}})
in mapfoldl at reduce.jl:357
in in at reduce.jl:402
in evaluate at /Users/goretkin/.julia/v0.3/ApproxFun/src/Spaces/Modifier/VectorSpace.jl:110
in getindex at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/Fun.jl:76
Another anomaly which the user may find difficult:
x = Fun(identity)
f = exp(x)./(1-x.^2)
f = exp(x)./(1-x.^2).^1
f = exp(x)./(1-x.^2).^1.0
I am following the README. The command countour(u)
doesn't work, as countour is undefined.
I assume I have to load a graphics package? Which one? I tried PyPlot, but this outputs an error starting with
ERROR: PyError (PyObject_Call) <type 'exceptions.TypeError'>
TypeError('float() argument must be a string or a number',)
I am using Julia 0.4-dev.
I just did Givens by hand, which slows down too much for large bandwidth. I'm sure this could be sped up by a large factor, either by using LA Pack on sub blocks, vectorization, or other tricks
The following no longer works, probably due to change in promotion
d=Disk()
u0 = Fun((x,y)->exp(-50x^2-40(y-.1)^2)+.5exp(-30(x+.5)^2-40(y+.2)^2),d)
B= [dirichlet(d) ,neumann(d)]
L=-lap(d)^2
h = 0.0005
timeevolution(2,B,L,u0,h)
Getting this to work would be helpful:
f(x,y) = exp(im*(x-y))
g = TensorFun(f)
h = LowRankFun(f)
f(0.5,0.5),g[0.5,0.5],h[0.5,0.5]
I'm trying to construct a vector-valued Fun
as follows, but it doesn't seem to be a proper Fun
. (It seems to work with FFun
, though)
julia> y0 = Fun(x->[1.0, 0.0])
Fun{VectorDomainSpace{2,UltrasphericalSpace{0},Float64},Float64}([1.0],VectorDomainSpace{2,UltrasphericalSpace{0},Float64}(UltrasphericalSpace{0}(Interval{Float64}(-1.0,1.0))))
julia> y0[1.0]
ERROR: BoundsError()
in clenshaw at /Users/goretkin/.julia/v0.3/ApproxFun/src/LinearAlgebra/clenshaw.jl:28
in evaluate at /Users/goretkin/.julia/v0.3/ApproxFun/src/Multivariate/VectorFun.jl:69
in getindex at /Users/goretkin/.julia/v0.3/ApproxFun/src/Fun/Fun.jl:69
Null(Operator) is slow and needs to be optimized. The culprit appear to be
3293 ...c/Operators/null.jl; applygivens!; line: 10
Right now norm first calculates sum(f^2) then takes sqrt. This will sometimes report 0.0, which seems misleading.
This is more of a goal than an issue. But ultimately there is no reason I can think of that the user shouldn't be able to use ApproxFun in BigFloat in the future (It would be BigFun!). Fun constructors, adaptive QR, and most other routines are based on a tolerance which is in turn based on an epsilon, usually eps(). If it were set to eps(T), then the routines would perform accordingly. Since the Chebyshev transform is based on the FFT, and Julia's FFT routines come from fftw, a step toward the BigFloat Chebyshev transform is included #53.
When I try a highly oscillatory function the IFun constructor oversamples, leaving a large tail of coefficients around machine precision. For instance,
f = IFun(x->sin(400*pi*x),[-1,1]);
println(length(f.coefficients))
8193
The actual length of the Chebyshev expansion should be a little larger than 400*pi, which is about 1250. The issue is probably in zerocfsIFun
(I think) because the resolution tolerance does not depend on the function. Somehow, the tolerance needs to depend on the absolute condition number of evaluating the function. See, for instance, here for the relative condition number:
http://en.wikipedia.org/wiki/Condition_number#One_variable
A quick estimate of the function's derivative (by a coarse finite difference stencil) suffices in the Chebfun constructor.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.