by

## Where communities thrive

• Join over 1.5M+ people
• Join over 100K+ communities
• Free without limits
##### Activity
Sheehan Olver
@dlfivefifty
I don't think so: that needs to be refactored. FullPDETest.jl suffered a bit of version control neglect.
On Tuesday let's discuss what an FEM analogue would be. That's a very different viewpoint from ApproxFun as one works with weak formulations, not operators.
Stefanos Carlström
@jagot
Yes, I realized I have too many questions to work efficiently at the moment
Sheehan Olver
@dlfivefifty
Domains.jl originated with @daanhb so some of the design is different. It also uses the term "Space" as things like R are Euclidean spaces. The analogue of AnyDomain() would be FullSpace{Any}().
Some more thought should be put in to how to unify this notion of Space with ApproxFuns.
Stefanos Carlström
@jagot
I see. I was also wondering, how I, if I knew a set of points x,y, would use transform to get the expansion coefficents, i.e. the best possible approximation/interpolation.
Some notion on how to perform eigendecompositions of operators acting on a certain basis, should also be encoded. For e.g. BSplines, since the basis functions are not orthogonal, you have an overlap matrix and hence you have a generalized eigenvalue problem.
Sheehan Olver
@dlfivefifty
Stefanos Carlström
@jagot
Maybe I do, but it seems that points then need to be attached to an existing Space.
Ok, Vandermonde seems to be what I want
Kirill Ignatiev
@ikirill
Is this supposed to work?
julia> dot(Fun(), Fun())
ERROR: StackOverflowError:
Stacktrace:
[1] dot(::Function, ::Function) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/generic.jl:665 (repeats 39122 times)
[2] top-level scope at none:0
Sheehan Olver
@dlfivefifty
Can you file an issue?
Kirill Ignatiev
@ikirill
Am I going crazy, or is there something wrong with transposes? JuliaApproximation/ApproxFun.jl#624
Up there when I was solving the associated Legendre ODE:
is there a way to get the diagonal operator whose entries are 1/sqrt(dot(psi[j], psi[k])), psi being the basis of a space?
It would make the self-adjoint ODE actually have a symmetric matrix when solving the generalized eigenvalue problem, but I couldn't find this.
Sheehan Olver
@dlfivefifty
It's just an operator that hasn't been added yet. Probably there should be a Normalized(Jacobi(1,1)) space and this done as a conversion operator.
(Before you ask "why not always use normalized", its because the unnormalized versions have rational formulae which are much faster.)
Mikael Slevinsky
@MikaelSlevinsky
I showed how to get a symmetric-definite + banded GEP for associated Legendre functions here (section 3.1 https://arxiv.org/abs/1711.07866). I have some ApproxFun-like code that I wrote to verify the correctness of the formulas, (but not the FMM speedup)
Mikael Slevinsky
@MikaelSlevinsky
Sorry, that should say "one way" to get banded symmetry
if I have 2D cartesian grid data at a certain resolution, is there any way to make a high precision Chebychev interpolant using ApproxFun?
ah, I see a comment above that answers this!
Soham Mukherjee
@soham1112
Hi, I was just trying out an example from https://github.com/JuliaApproximation/ApproxFunExamples/blob/master/PDEs/Rectangle%20PDEs.ipynb. When trying the example for solving the convection equation (see In[10]), I end up getting an error if I end up changing the interval for t, say, for example, from 0 .. 2 to 2 .. 3. ERROR: LoadError: LoadError: Cannot convert coefficients from Chebyshev(-1.0..1.0) to Chebyshev(-1.0..1.0)⊗ConstantSpace(Point(2.0))
I don't understand why this should be giving an error, since I should be able to choose an arbitrary time-interval. Am I missing something?
Sheehan Olver
@dlfivefifty
It’s using an embedding of a 1D interval in the 2D plane where the interval is assumed to be at t = 0. You can use a fully 2D version: u0 = Fun((x,_) -> ..., Chebyshev(-1.0..1.0)⊗ConstantSpace(Point(2.0)))
Soham Mukherjee
@soham1112
Ah! I see. Thanks--that worked like a charm! However, I'm curious if there's a reason for this default behaviour?
Mitko Georgiev
@mitkoge
Hi,
I tried tzhe following
using ApproxFun
sp1= Chebyshev(0.0..50.0) ∪ Chebyshev(50.0..60.0)
x1=Fun(sp1); g1 = Fun(exp(-x1^2),sp1)
s1= sum(g1) #NaN
sp2= Chebyshev(0.0..5.0) ∪ Chebyshev(5.0..60.0)
x2=Fun(sp2); g2 = Fun(exp(-x2^2),sp2)
s2= sum(g2) #0.8862269254527582
i expected s1~= s2. but it seems it is not.
Mikael Slevinsky
@MikaelSlevinsky
I pushed a fix for this on this branch JuliaApproximation/ApproxFun.jl@07f50d0. If the tests pass, then it could be merged.
Mitko Georgiev
@mitkoge
Wow i am impressed by your reactivity! Fixed, pushed and merged in a day! Thank you!
Mitko Georgiev
@mitkoge

I try to reproduce an example from this chat given by @cortner

using ApproxFun
dom = Interval(0.001, 1) * PeriodicInterval(-pi, pi)
space = Space(dom)
Dr = Derivative(space, [1,0])
Dθ = Derivative(space, [0,1])
Mr = Multiplication(Fun( (r, θ) -> r, space ), space)
rDr = Mr * Dr
L = rDr * rDr + Dθ * Dθ

but PeriodicInterval(-pi, pi) is no more recognized.
What alternative can be used? Circle?

Sheehan Olver
@dlfivefifty
It’s been renamed PeriodicSegment as it’s really an oriented line segment
Mitko Georgiev
@mitkoge

Thank you,
i tried

julia> dom = Interval(0.001, 1) * PeriodicSegment(-pi, pi)
ERROR: MethodError: no method matching *(::Interval{:closed,:closed,Float64}, ::PeriodicSegment{Floa
t64})

and

julia> dom = Segment(0.001, 1) * PeriodicSegment(-pi, pi)
ERROR: MethodError: no method matching *(::Segment{Float64}, ::PeriodicSegment{Float64})
Sheehan Olver
@dlfivefifty
Ah and * has changed to imes
imes
\times
Mitko Georgiev
@mitkoge

thank you.

using ApproxFun,LinearAlgebra
dm = Segment(0.001,1) × PeriodicSegment(-pi,pi)
sp = Space(dm)
Dr = Derivative(sp, [1,0]); Dθ = Derivative(sp, [0,1])
Mr = Multiplication(Fun( (r, θ) -> r, sp ), sp)
rDr = Mr * Dr
Lr = rDr * rDr; L = Lr + Dθ * Dθ

i modified it and now can play with.

Simon Etter
@ettersi

Not sure whether this is the right place to ask since technically my question concerns SingularIntegralEquations.jl, but that package seems to be related enough to ApproxFun.jl that I'll go ahead anyway.

In the code of SingularIntegralEquations it says

# stieltjesintegral is an indefinite integral of stieltjes
# normalized so that there is no constant term

Can someone elaborate on that normalisation condition?

Sheehan Olver
@dlfivefifty
Sure, fine to ask here. This normalisation condition comes from the fact that Jacobi polynomials have nice integration formulae coming from the weighted differentiation d/dx[(1-x)^a *(1+x)^b P_n^(a,b)] = C (1-x)^{a-1} (1+x)^{b-1} P_{n+1}^(a-1,b-1). Roughly speaking, for functions that vanish at ±1, derivatives and Stieltjes transforms commute, so provided a,b > 0 we can use this to write the integral of a Stieltjes transforms in terms of the integral of the integrand.
Note that the n = 0 case is special, and dictates the normalisation constant (as the construction for other n will automatically decay something like O(z^(-n)) )
For the cases actually implemented, I think its normalized to go to zero at ∞.
Sorry, that's not right.
You always have logarithmic growth (unless the integrand integrates to zero).
Sheehan Olver
@dlfivefifty
On -1..1 it just vanishes apart from the log term:
julia> x = Fun();

julia> f = exp(x)/sqrt(1-x^2);

julia> z = 10_000; stieltjesintegral(f,z) - log(z) * sum(f)
-0.00017756097918919522
julia> f = exp(x) * sqrt(1-x^2);

julia> z = 10_000; stieltjesintegral(f,z) - log(z) * sum(f)
-4.264886882765495e-5
Oh, wait, that's true for other intervals too:

julia> x = Fun(1..2);

julia> f = exp(x) * sqrt((x-1)*(2-x));

julia> z = 100_000; stieltjesintegral(f,z) - log(z) * sum(f)
-2.8356239955229512e-5
Sorry for making it seem more complicated than it actually it is. It's really simple: it's normalized so that stieltjes(f,z) = log(z) *sum(f) + o(1).