JuliaApproximation/ApproxFun.jl

Sheehan Olver

@dlfivefifty

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).

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

Ashton Bradley

@AshtonSBradley

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.

Ashton Bradley

@AshtonSBradley

ok thanks, good to know

Ashton Bradley

@AshtonSBradley

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?

Ashton Bradley

@AshtonSBradley

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?

Where communities thrive

People

Repo info

Activity