using Plots; pyplot()
using Blink, Interact
using ApproxFun
plot(rand(3))
function chebpts(N)
return cos.((0:N)*pi/N)
end
chebpts(5)
"""
Compute and return the Chebyshev differentiation matrix.
Based on Trefethen, Spectral Methods in Matlab, page 54, cheb.m
function [D,x] = cheb(N)
if N==O, D=O; x=1; return, end
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1); 2].*(-1)."(O:N)'; X = repmat(x,1,N+1);
dX = X-X';
D = (c*(1./c)')./(dX+(eye(N+1)));
D = D - diag(sum(D'));
"""
using LinearAlgebra
function chebdiff(N)
if N==0
D=zeros(1,1)
x=ones(1.)
else
x = chebpts(N)
c = [2; ones(N-1); 2] .* (-1).^(0:N)
X = repeat(x,outer=[1,N+1])
dX = X - X'
D = (c*(1 ./ c)')./(dX + I)
D = D - diagm(vec(sum(D,dims=2)))
end
return D,x
end
chebdiff(5)
chebdiff(2)[1]
w = Window()
f(x) = exp.(sin.(5 .* x)).*cos.(x)
af = Fun(f,-1..1)
ui = @manipulate for N=20:50
x = chebpts(N)
D = chebdiff(N)[1]
v = f(x)
plot(af,linewidth=5.,color=colorant"lightgrey", label="f")
plot!(af',linewidth=5.,color=colorant"lightgrey", label="f'")
plot!(x,v,marker=:x,label="f at chebpts")
plot!(x,D*v,marker=:o,label="Dv at chebpts")
end
body!(w,ui)
a = -5;
y0 = 1.
N = 11
D,x = chebdiff(N)
u = zeros(N+1)
D1 = D[1:end-1,1:end-1]
d0 = D[1:end-1,end]
u1 = (D1-a*I) \(-y0*d0)
u[end] = y0
u[1:end-1] = u1
plot(Fun(x -> exp(a*(x+1))),linewidth=5.,color=colorant"lightgrey")
plot!(x,u)
A = 25*[-0.01 1; -1 -0.01]
y0 = [0.;1.]
N = 50
x = chebpts(N)
D = chebdiff(N)[1]
U = zeros(N+1,2)
D1 = D[1:end-1,1:end-1]
d0 = D[1:end-1,end]
u2 = (kron(I(2),D1)-kron(A,I(N)))\(-kron(y0,d0))
U[1:end-1,:] = u2
U[end,:] = y0
plot(x,U)
using Printf
ui = @manipulate for N=10:250
A = 50*[-0.01 1; -1 -0.01]
y0 = [0.;1.]
x = chebpts(N)
D = chebdiff(N)[1]
U = zeros(N+1,2)
D1 = D[1:end-1,1:end-1]
d0 = D[1:end-1,end]
u2 = (kron(I(2),D1)-kron(A,I(N)))\(-kron(y0,d0))
U[1:end-1,1:end] = u2
U[end,:] = y0
plot(x,U)
plot!(title=@sprintf("N=%i",N))
end
body!(w,ui)
A = 25*[-0.01 1; -1 -0.01]
y0 = [0.;1.]
N = 200
x = chebpts(N)
D = chebdiff(N)[1]
U = zeros(N+1,2)
D1 = D[1:end-1,1:end-1]
d0 = D[1:end-1,end]
u2 = (kron(I(2),D1)-kron(A,I(N)))\(-kron(y0,d0))
U[1:end-1,:] = u2
U[end,:] = y0
plot(x,U)
#[U[1,:]; (expm((A/25)*50)*y0)' ]
using LinearAlgebra
N = 300
x = chebpts(N)
Ue = zeros(N+1,2)
A = [-0.01 1; -1 -0.01];
for i=1:N
Ue[i,:] = (exp(A*(x[i]+1)*25)*y0)'
end
plot(x,sum(Ue.^2,dims=2))
A