bilaplacian

[rveltz]

Hi,

I have trouble implementing a bilaplacian with Neumann boundary condition. First, for the linear operator

d=Interval()^2 op = lap(d)^2

I get:

ERROR: AssertionError: b >= a in call at /Users/rveltz/.julia/v0.4/ApproxFun/src/Spaces/Ultraspherical/UltrasphericalOperators.jl:237 in promotedomainspace at /Users/rveltz/.julia/v0.4/ApproxFun/src/Operators/spacepromotion.jl:111 in promotedomainspace at /Users/rveltz/.julia/v0.4/ApproxFun/src/PDE/KroneckerOperator.jl:105 in promotedomainspace at /Users/rveltz/.julia/v0.4/ApproxFun/src/Operators/banded/algebra.jl:578 in promotedomainspace at /Users/rveltz/.julia/v0.4/ApproxFun/src/Operators/spacepromotion.jl:106 in promotetimes at /Users/rveltz/.julia/v0.4/ApproxFun/src/Operators/banded/algebra.jl:352 in * at /Users/rveltz/.julia/v0.4/ApproxFun/src/Operators/banded/algebra.jl:459 in power_by_squaring at intfuncs.jl:78 in ^ at intfuncs.jl:108

Can you show me how to define lap(d)^2 please?

Thank you for your help.

[dlfivefifty]

PDE solving is very experimental at the moment. This example can be made to work by the following:

Pkg.checkout(“ApproxFun”,”development”)     # get new bug fixes

dx=dy=Interval() 
d=dx*dy 
Dx=Derivative(dx);Dy=Derivative(dy) 
L=Dx^4⊗I+Dx^2⊗Dy^2+I⊗Dy^4                 # equivalent to lap(d)^2

 K=kronfact([dirichlet(d); neumann(d); L],100,100)   # construct 100x100 matrix

# create a RHS
 x=Fun(identity,dx);y=Fun(identity,dy) 
 G=[real(exp(-1+1.im*y)); real(exp(1+1im*y)); real(exp(x-1im)); real(exp(x+1im))] 

# solve for u
u=K\G 

 using Plots 
surface(u)                  # plot the solution in PyPlot

On 22 Feb 2016, at 10:41 PM, rveltz notifications@github.com wrote:

Hi,

I have trouble implementing a bilaplacian with Neumann boundary condition. First, for the linear operator

d=Interval()^2 op = lap(d)^2

I get:

ERROR: AssertionError: b >= a in call at /Users/rveltz/.julia/v0.4/ApproxFun/src/Spaces/Ultraspherical/UltrasphericalOperators.jl:237 in promotedomainspace at /Users/rveltz/.julia/v0.4/ApproxFun/src/Operators/spacepromotion.jl:111 in promotedomainspace at /Users/rveltz/.julia/v0.4/ApproxFun/src/PDE/KroneckerOperator.jl:105 in promotedomainspace at /Users/rveltz/.julia/v0.4/ApproxFun/src/Operators/banded/algebra.jl:578 in promotedomainspace at /Users/rveltz/.julia/v0.4/ApproxFun/src/Operators/spacepromotion.jl:106 in promotetimes at /Users/rveltz/.julia/v0.4/ApproxFun/src/Operators/banded/algebra.jl:352 in * at /Users/rveltz/.julia/v0.4/ApproxFun/src/Operators/banded/algebra.jl:459 in power_by_squaring at intfuncs.jl:78 in ^ at intfuncs.jl:108

Can you show me how to define lap(d)^2 please?

Thank you for your help.

— Reply to this email directly or view it on GitHub https://github.com/ApproxFun/ApproxFun.jl/issues/298.

[rveltz]

Thank you for your answer, it works very well.

May I ask for another related hint? How can you specify different boundary conditions (than dirichlet(d)), for example a linear function of u and Du on the boundary?

Thank you,

Best.

[dlfivefifty]

Eventually you’ll be able to do [dirichlet(d)+neumann(d);dirichlet(d)-neumann(d)], for example, but that doesn’t work at the moment. For now you can write it out explicitly in Kronecker form:

K=kronfact([(ldirichlet(dx)+lneumann(dx))⊗I; (rdirichlet(dx)+rneumann(dx))⊗I; I⊗(ldirichlet(dy)+lneumann(dy)); I⊗(rdirichlet(dy)+rneumann(dy)); (ldirichlet(dx)-lneumann(dx))⊗I; (rdirichlet(dx)-rneumann(dx))⊗I; I⊗(ldirichlet(dy)-lneumann(dy)); I⊗(rdirichlet(dy)-rneumann(dy)); L],100,100)

u=K\G

On 23 Feb 2016, at 4:39 AM, rveltz notifications@github.com wrote:

Thank you for your answer, it works very well.

May I ask for another related hint? How can you specify different boundary conditions (than dirichlet(d)), for example a linear function of u and Du on the boundary?

Thank you,

Best.

— Reply to this email directly or view it on GitHub https://github.com/ApproxFun/ApproxFun.jl/issues/298#issuecomment-187283802.

[dlfivefifty]

PS There's a major revamp of PDEs in the works, so this won't be simplified until after that is done. But glad to hear the package is useful!

[rveltz]

Hi,

I am coming back to you because I am not sure I understand if your answer apply only to PDE. Hence, is it difficult to implement a boundary condition like u(0) = B(u), where u solves a 1d integro-differential problem L? (This question is also related to my previous questions about Volterra equations).

[dlfivefifty]

I'm not sure I understand, but 1D is much more developed than PDEs. Your example should eventually work, if you give us a concrete example we can have a go

Sent from my iPhone

On 9 Mar 2016, at 18:20, rveltz notifications@github.com wrote:

Hi,

I am coming back to you because I am not sure I understand if your answer apply only to PDE. Hence, is it difficult to implement a boundary condition like u(0) = B(u), where u solves a 1d integro-differential problem L? (This question is also related to my previous questions about Volterra equations).

— Reply to this email directly or view it on GitHub.

[rveltz]

For example solve for x>=0 0 = -u(x)' - f(x)u(x) with the constraint u(0) = integral (Phi(x) * u(x), x=0..infinity)

[dlfivefifty]

What’s Phi and f?

On 9 Mar 2016, at 10:32 PM, rveltz notifications@github.com wrote:

For example solve for x>=0 0 = -u(x)' - f(x)u(x) with the constraint u(0) = integral (Phi(x) * u(x), x=0..infinity)

— Reply to this email directly or view it on GitHub https://github.com/ApproxFun/ApproxFun.jl/issues/298#issuecomment-194252485.

[rveltz]

two functions that are not L1. f=x->x^2 and Phi=x->0.4*x+f(x)

[dlfivefifty]

Isn’t the solution u = 0?

On 10 Mar 2016, at 6:56 PM, rveltz notifications@github.com wrote:

two functions that are not L1. f=x->x^2 and Phi=x->0.4*x+f(x)

— Reply to this email directly or view it on GitHub https://github.com/ApproxFun/ApproxFun.jl/issues/298#issuecomment-194721198.

[dlfivefifty]

In any case, the following works (or you can also truncate the domain and use Chebyshev):

S=JacobiWeight(0.,4,Ray())
D=Derivative(S)
x=Fun(domain(S))
f=x^2
Φ=0.4x+f

B=DefiniteIntegral()[Φ]
u=[B;
 D+f]\[1.]

[rveltz]

Hi,

Sweeet!

Sure, u=0 is solution ; I want the eigenvalues of the system.

• if I truncate, say S=Interval(0,1), '[B;D+f]' gives u(0) = Bu and u(1)=Bu? can I select only u(0) = Bu?
• if I pass your code to eigvals, it will look for eigenvalues lambda solving lambda u = u'+fu but for which boundary condition? u(0) = Bu or (lambda-1)u(0) = Bu ?

Thank you for your response, I am very close now to my needs.

[dlfivefifty]

It does u(0), not (λ-1)u(0).

I added eigs, which now prunes out false eigenvalues by looking at the tail of the eigenvectors:

S=Chebyshev()
A=[dirichlet(S);Derivative(S,2)]
λ,V=eigs(A,50)

Here are your examples, both on the truncated interval

S=Chebyshev([0.,10.])
D=Derivative(S)
x=Fun(domain(S))
f=x^2
Φ=0.4x+f

B=DefiniteIntegral()[Φ]
λ,V=eigs([B;D+f],100)

plot(V[1])

or on the ray

S=JacobiWeight(0.,4,Ray())
D=Derivative(S)
x=Fun(domain(S))
f=x^2
Φ=0.4x+f
B=DefiniteIntegral()[Φ]
λ,V=eigs([B;D+f],200)

plot(0.:.01:10.,V[1])

[rveltz]

That's great!! Thank you a lot.

Something is still not clear about the boundary conditions. When you write λ,V=eigs([B;D+f],200) for the case S=Chebyshev([0.,10.]) does it mean the BCs u(0) = Bu and u(10)=Bu or just one of them?

How can I select say the left one like when doing ldririchlet(S)?

[dlfivefifty]

Here it's just doing

\int_0^10 u(x) Φ(x) dx = 0

If you want

\int_0^10 u(x) Φ(x) dx = u(0)

you should define

B=DefiniteIntegral()[Φ]-ldirichlet()

[rveltz]

Ah! I got it now, this was not very clear.

Thank you a lot for your help,

[dlfivefifty]

Sorry it wasn't clear. Now that you understand how it's meant to work, please feel free to contribute to the Wiki documentation.

I'm reopening, since the original issue about bilaplacian is not fixed yet.

[rveltz]

Hi,

To come back on your function eigs, I am wondering if you could feed your evaluation function of the linear operator to a krylov eigensolver. If I understand what you are doing, you 'just' compute the matrix of the linear operator in the basis of the first n polynomials and give it to eigs. But this might not be well adapted to the eigenproblem, e.g. whether you want the rightmost ev for example.

This would led to a more precise computation of the spectrum,

Best regards,

[dlfivefifty]

Yes the current implementation is a quick hack that does 'just' what you said.

Down the line we might do simultaneous inverse iteration using \, which would have better adaptive properties.

I think Krylov methods are better off leaving to the user.

Sent from my iPhone

On 15 Mar 2016, at 04:10, rveltz notifications@github.com wrote:

Hi,

To come back on your function eigs, I am wondering if you could feed your evaluation function of the linear operator to a krylov eigensolver. If I understand what you are doing, you 'just' compute the matrix of the linear operator in the basis of the first n polynomials and give it to eigs. But this might not be well adapted to the eigenproblem, e.g. whether you want the rightmost ev for example.

This would led to a more precise computation of the spectrum,

Best regards,

— Reply to this email directly or view it on GitHub.

[dlfivefifty]

The following works (soon to be merged):

dx=dy=Interval()
d=dx*dy
Dx=Derivative(dx);Dy=Derivative(dy)
L=Dx^4⊗I+2*Dx^2⊗Dy^2+I⊗Dy^4

A=[dirichlet(dx)⊗eye(dy);
        eye(dx)⊗dirichlet(dy);
        neumann(dx)⊗eye(dy);
        eye(dx)⊗neumann(dy);
         L]

F=[Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A[1]));
    Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A[2]));
    Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A[3]));
    Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A[4]));
    Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A[5]));
    Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A[6]));
    Fun((x,y)->-imag(exp(x+1.0im*y)),rangespace(A[7]));
    Fun((x,y)->-imag(exp(x+1.0im*y)),rangespace(A[8]));
    0]

u=linsolve(A,F;tolerance=1E-10)