by

## Where communities thrive

• Join over 1.5M+ people
• Join over 100K+ communities
• Free without limits
##### Activity
Sheehan Olver
@dlfivefifty
Which happens to be consistent with logkernel (whose normalization is imposed by the actual definition):
julia> logkernel(f, z) -log(abs(z))/π * sum(f)
-9.026071512430178e-6
Simon Etter
@ettersi
If I understand correctly, I can summarise what you said as: stieltjes(f,z) = integral log(z - w) f(w) dw + C with Csuch that stieltjes(f,z) = log(z) *sum(f) + o(1) for z to Inf. Correct?
Sheehan Olver
@dlfivefifty
Yes
Simon Etter
@ettersi
Great, thanks a lot!
Mitko Georgiev
@mitkoge

Hi i was playing with modified example from Quantum states.jl

Using  ApproxFun
fc1 = Fun(1, 0..1); fc2 = Fun(2, 1..2); fc3 = Fun(3, 2..3)

fc= fc1+fc2+fc3
D = Derivative(); sp = space(fc); x= Fun(sp);
B = [Dirichlet(sp); continuity(sp,0:1)]
λ, v = ApproxFun.eigs(B, -D^2 + fc, 500,tolerance=1E-10)
#returns result

λ, v = ApproxFun.eigs(B, -D^2 - 1/x*D + fc, 500,tolerance=1E-10)
# says Not implemented  but if fc=fc1+fc2 it manages to return result

Is 1/x*D kind of a problematic composition?

Sheehan Olver
@dlfivefifty
Yes. It would be easy to work around by multiplying by x and updating the generalised eigenvalue call appropriately in the implementation of eigs.
Mitko Georgiev
@mitkoge

Thanks,
you mean calling like this

λ, v = ApproxFun.eigs_overx(B, -x*D^2 - D + x*fc, 500,tolerance=1E-10)

and somehow divide C1 to x before calling to

    λ,V = eigen(A1,C1)

inside the copy of eigs?

Sheehan Olver
@dlfivefifty
No multiply by x: instead of solving A*v = λ*C*v we solve x * A * v = λ * x * C* v. That is, in line 34 change to C = Conversion(ds,rangespace(A)); C = x * C should be enough.
Note eigs is quite stale and needs to be rewritten, but that's low priority for me. I'll accept any PRs.
Mitko Georgiev
@mitkoge
Thanks for the great support,
Mitko Georgiev
@mitkoge

i followed you and copy paste defined

    function eigs(Bcs_in::Operator,A_in::Operator,fmul::Fun,n::Integer;tolerance::Float64=100eps())

with the correction

        C = Conversion(ds,rangespace(A)); C = fmul * C

now have only to check if it works correct. :-)

λ, v = ApproxFun.eigs(B, L, x, 500,tolerance=1E-10)
Mitko Georgiev
@mitkoge
:point_up: 22. November 2018 12:35
I intendet to help with replacing it in the docs but as i see it is already fixed.
Unfortunately it might be not at http://juliaapproximation.github.io/ApproxFun.jl/latest
is it possible to solve a nonlinear PDE in approxfun? for example I would like to add a simple nonlinear term to this problem
https://github.com/JuliaApproximation/ApproxFun.jl#solving-partial-differential-equations
, presumably switching to newton solver, but not sure how to inform it about boundary conditions etc in 2D
Sheehan Olver
@dlfivefifty
Not yet: I want to make the PDE solver faster before trying to implement that sort of thing. For now you can try doing Newton iteration by hand.
ok thanks, good to know
for the basic procedure of setting up for manual iteration, can I take it that this https://github.com/JuliaApproximation/ApproxFunExamples/blob/master/ODEs/Blasius%20and%20Falkner-Skan.jl
is a reasonable starting point, that I just need to generalise to 2D?
Ok, maybe ignore that
Sheehan Olver
@dlfivefifty
It’s pretty simple: just differentiate wrt u but treating d/dx and d/dy as constants.
To get the Jacobian
Mitko Georgiev
@mitkoge
excuse if my question is too naive but i try this
x= Fun(0..1); sx=space(x);
a=Fun(0..pi); sa=space(a);
fa= DefiniteIntegral(x*a,sx)
# i expect fa= a here
but i get MethodError. Are integrals of two arguments possible?
Mitko Georgiev
@mitkoge
or may be like this
x= Fun(0..1); sx=space(x)
a=Fun(0..pi); sa=space(a)
fm= x*a; sm= space(fm)
fI=  DefiniteIntegral(sx)
(fI*fm)(0.2)
#expect 0.2
Sheehan Olver
@dlfivefifty

Are you expecting it to be a 2D function? And do you actually expect the result to be 0.1 (that is, a/2 since you integrate out x)?

That's possible via:

x,a = Fun((0..1) × (0..π))
Q = DefiniteIntegral(0..1) ⊗ I
(Q*(x*a))(0.0,0.2)
Mitko Georgiev
@mitkoge
Wow exciting how elegant is this! Thank you!
As you can see i am exploring it the dummy way so please bear with me.
Mitko Georgiev
@mitkoge
let Q2= Q (x a) . Q2 is still 2D Fun((x,a)-> smth). Is there a way to curry it to 1D function Q1 that is Fun(a-> smth)?
Sheehan Olver
@dlfivefifty
Not an easy way. What we actually want is DefiniteIntegral(S, [1,0]) that wraps the kronecker product. It’s not too hard to add this but I don’t have time right now
Mitko Georgiev
@mitkoge

to practice your guiding i wrote

function Hankel(fx::Fun, spx,spa::Space)
spxa= spx ⊗ spa
fx2= Fun((x,a)-> fx(x), spxa)
fj= Fun((x,a) -> besselj0(x*a), spxa)
fa2= (DefiniteIntegral(dmx1) ⊗ I)*(fx2*fj*x)
fa= Fun(a->fa2(0,a), spa)
end

from the definition. It seems to work. Thank you!
Is this the style ApproxFun is expected to be used?

Sheehan Olver
@dlfivefifty

Not an easy way.

I take that back. I think you can always do ApproxFun.SpaceOperator(op, newdomainspace, newrangespace) change the spaces of an operator. So just write Q with the "right" spaces

Is this the style ApproxFun is expected to be used?

In Julia, capital letters are only used when it returns a special type of the same name. So this would be better as lower case. Also, there is support for Green's functions in https://github.com/JuliaApproximation/SingularIntegralEquations.jl including Helmholtz / hankel kernels

Mitko Georgiev
@mitkoge
I see. I will check SpaceOperator and will have to learn what rangespace is.
And thank you for the link. It is probably some more efficient implementation? Will check it.
I tried hankel as an learning example.
The kind of high level style the ApproxFun works is fascinating me.
Like a poetry for the crowd. I do not necessarily understand but like it. :-)
Mitko Georgiev
@mitkoge
here is the revised
function hankel(fx::Fun, spx,spa::Space)
spxa= spx ⊗ spa
fx2= Fun((x,a)-> fx(x), spxa)
fjx= Fun((x,a)-> x*besselj0(x*a), spxa)
fa2= (DefiniteIntegral(spx.domain) ⊗ I)*(fx2*fjx)
fa= Fun(a->fa2(0,a), spa)
end
Mitko Georgiev
@mitkoge

when i blind try

spx= Space(0..1); spa= Space(0..π); spxa= spx ⊗ spa
Q2= DefiniteIntegral(spx.domain) ⊗ I
Q1= ApproxFun.SpaceOperator(Q2, spxa, spa )

x=Fun((x,a)->x,spxa); a=Fun((x,a)->a,spxa)
(Q2*(x*a))(0,1)  #=0.5
(Q1*(x*a))(1)      #=0.5*pi/2

it seems to work up to scaling cofficient

Mitko Georgiev
@mitkoge
i gues it is because the default domain for I is -1..1
Mitko Georgiev
@mitkoge
may i ask if a least quare fit by bivariate Fun to data is possible?
Given data is Array{Float64,2}.
I liked your nice 1D example from the docs and wonder how to extend to 2D.
Sheehan Olver
@dlfivefifty
Yes an example is in the FAQ
Mitko Georgiev
@mitkoge
Thank you! I found it and will have a look.
Is the example for standard grid ponts also extendable for 2D?
Or it is not preferable because data points in 2D are mostly large ammont?
Art Gower
@arturgower

I'm quite sure Fun used to work on intervals in the complex domain. But now

f(x) = cos(x)
Fun(f, Interval(1.0+1.0im, 2.0+2.0im))

throws the error

ERROR: MethodError: no method matching isless(::Complex{Float64}, ::Complex{Float64})
Closest candidates are:
isless(::Missing, ::Any) at missing.jl:66
isless(::InfiniteArrays.OrientedInfinity{Bool}, ::Number) at /.julia/packages/InfiniteArrays/Z4yap/src/Infinity.jl:145
isless(::Number, ::InfiniteArrays.OrientedInfinity{Bool}) at /.julia/packages/InfiniteArrays/Z4yap/src/Infinity.jl:144
...
Stacktrace:
[1] <(::Complex{Float64}, ::Complex{Float64}) at ./operators.jl:260
[2] >(::Complex{Float64}, ::Complex{Float64}) at ./operators.jl:286
[3] isempty(::Interval{:closed,:closed,Complex{Float64}}) at /.julia/packages/IntervalSets/xr34V/src/IntervalSets.jl:153
Sheehan Olver
@dlfivefifty
Use Segment(a,b) for line segments in the complex plane. (The previous support for intervals in the complex plane violated the definition of an interval.)
Art Gower
@arturgower
A ha! Yes I see. Thanks so much
Shi Pengcheng
@shipengcheng1230

Hello, I was playing around with the poisson equation example and wonder if I could replace the RHS $f$ with something like $\delta (x) \delta (y)$. I tried to construct the RHS like this:

fx = KroneckerDelta()
fy = KroneckerDelta()
f = Fun((x,y) -> fx(x) * fy(y))

But I got the error that ERROR: MethodError: no method matching isless(::Int64, ::Nothing). What is the proper way for me to do that? Thanks in advance!

Sheehan Olver
@dlfivefifty
Essentially you want to calculate the Greens function? That’s a tough question which we are looking at right now.
Shi Pengcheng
@shipengcheng1230
Thanks for the reply. I guess I have to wait for now :)
Quentin C.
@robocop
Hello ! I would like to consider ApproxFun for the following problem: I have a N-L operator Phi: L_s -> L_s, where L_s is the set of continuous functions that goes exponentially fast to zero at speed s > 0, that is: x belongs to Ls if \sup{t \geq 0} |x(t)| e^{s t} < +oo. Moreover given x in L_s, Phi(x) is the solution of a Volterra integral equation of the form Phi(x)(t) = K_0(x)(t) + \int_0^t{K(x) (t, u)Phi(x) (u) du}, where the Kernels K_0(x) and K(x) are known explicitly in term of x. I would like to compute the spectrum of the Frechet differential of Phi at the point x = 0. Do you think it is possible to do that?
Let me know if the question is unclear... I do not have a background in numerical simulations. Thanks :)
Quentin C.
@robocop
I think I have started to understand how to do that with ApproxFun. But, I have got few questions. How to encode the domain $\mathbb{R}_+$?
and how to encode $\{ (x, y) \in mathbb{R}_+, x \geq y \}$?
Quentin C.
@robocop
More generally the fact that I want to work on $\mathbb{R}_+$ and not on a segment seems to be a problem. Is that correct?
Sheehan Olver
@dlfivefifty
Sorry forgot to reply. I have a student working on Volterra integral equations so once that’s up and running your problem should be easy. We’d love to know the application so we can mention it in the paper
Quentin C.
@robocop
Ok great! Looking forward to see that. The Volterra equation I am looking at appears in a neuroscience problem (see https://arxiv.org/abs/1810.08562) The dynamic of the mean-field network can be reduce to a N-L Volterra equation (see equations (7) and (8) of the paper). Btw the same problem can be solved numerically by looking at the associated PDE (see equation (3)) - but I am curious to know if more efficient numerically methods can be develop by specifically exploit this Volterra equation.