Solving PDE's on quasi-periodic domains

[bhaveshshrimali]

Hi @dlfivefifty

Firstly, thanks for ApproxFun. :) !!

I wanted to translate a Burger's equation solver from chebfun to ApproxFun. Could you point me to the equivalent of spinop in ApproxFun, if it exists. Or more generally the way to achieve the following:

function u = Burgers(init, tspan, s, visc)

S = spinop([0 1], tspan);

S.lin = @(u) + visc*diff(u,2);
S.nonlin = @(u) - 0.5*diff(u.^2);
S.init = init;
u = spin(S,s,1e-4,'plot','off'); 

Thanks,

[dlfivefifty]

What does spin do?

[bhaveshshrimali]

Spin is the stiff pde integrator object in chebfun. https://www.chebfun.org/docs/guide/guide19.html

The whole idea is to solve Burger's equation, for u(t, x), say on (0, 1] x [0, 1] where the initial condition is a gaussian, u(0, x) and the boundary conditions are periodic in x, namely u(t, 0) = u(t, 1) and u'(t, 0) = u'(t, 1)...

A slightly related query that I also have is as follows -- in chebfun apparently you can do

uu = chebfun(c, [0 1], 'trig', 'coeffs');   % c are the coefficients 
u = chebfun(@(t) uu(t - 0.5), [0 1], 'trig');

where c is the set of coefficients used in the trigonometric / fourier expansion of the function. Is it possible to do something similar using ApproxFun?

The code I am currently playing with is

using ApproxFun
using LinearAlgebra
function GRF(N, m, γ, τ, σ, type)
    if type == "periodic"
        my_const = pi
    end

    my_eigs = (sqrt(2) * (abs(σ) .* ((my_const .* (1:N)').^2 .+ τ^2).^(-γ/2)))'
    
    println(size(my_eigs))
    ξₐ = randn(N,1); # vector
    α = my_eigs .* ξₐ; # matrix
    ξᵦ = randn(N,1);
    β = my_eigs .* ξᵦ;
    a = α ./ 2;
    b = β ./ 2;
    c = [reverse(a, dims=1) - reverse(b, dims=1) .* 1im; m .+ 0*1im; a + b .* 1im];
    uu = Fun(c, Fourier(0..1));   # <FAILS HERE
    # u = Fun(t->uu(t - 0.5), Fourier(0..1));

end

which fails at uu = Fun(c, Fourier(0..1)).

Thanks for your help (and sorry for asking all of this at once).

[bhaveshshrimali]

For the latter, maybe this works (I should have taken a look at the docs more carefully).. alteast no errors this time

uu = Fun(Fourier(0..1), [c...])
u = Fun(t->uu(t - 0.5), Fourier(0..1));

with the complete function being...

function GRF(N, m, γ, τ, σ, type)
    if type == "periodic"
        my_const = pi
    end

    my_eigs = (sqrt(2) * (abs(σ) .* ((my_const .* (1:N)').^2 .+ τ^2).^(-γ/2)))'
    
    println(size(my_eigs))
    ξₐ = randn(N,1); # vector
    α = my_eigs .* ξₐ; # matrix
    ξᵦ = randn(N,1);
    β = my_eigs .* ξᵦ;
    a = α ./ 2;
    b = β ./ 2;
    c = [reverse(a, dims=1) - reverse(b, dims=1) .* 1im; m .+ 0*1im; a + b .* 1im];
    uu = Fun(Fourier(0..1), [c...])
    u = Fun(t->uu(t - 0.5), Fourier(0..1));

end