Solving PDE's on quasi-periodic domains
bhaveshshrimali opened this issue · 3 comments
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,
What does spin do?
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).
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