In [6]:
using Plots; pyplot()
using Interact
using ApproxFun
plot(rand(3))
Out[6]:
In [7]:
function chebpts(N)
    return cos((0:N)*pi/N)
end
chebpts(5)
Out[7]:
6-element Array{Float64,1}:
  1.0     
  0.809017
  0.309017
 -0.309017
 -0.809017
 -1.0     
In [8]:
"""
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'));
"""    
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 = repmat(x,1,N+1)
        dX = X - X'
        D = (c*(1./c)')./(dX + eye(N+1))
        D = D - diagm(vec(sum(D,2)))
    end
    return D,x
end
chebdiff(5)
Out[8]:
(
6x6 Array{Float64,2}:
  8.5       -10.4721     2.89443   -1.52786    1.10557   -0.5     
  2.61803    -1.17082   -2.0        0.894427  -0.618034   0.276393
 -0.723607    2.0       -0.17082   -1.61803    0.894427  -0.381966
  0.381966   -0.894427   1.61803    0.17082   -2.0        0.723607
 -0.276393    0.618034  -0.894427   2.0        1.17082   -2.61803 
  0.5        -1.10557    1.52786   -2.89443   10.4721    -8.5     ,

[1.0,0.8090169943749475,0.30901699437494745,-0.30901699437494734,-0.8090169943749473,-1.0])
In [9]:
f(x) = exp(sin(5*x)).*cos(x)
af = Fun(f,[-1,1])
@manipulate for N=20:50
    x = chebpts(N)
    D = chebdiff(N)[1]
    v = f(x)
    plot(af,linewidth=5.,color=colorant"lightgrey")
    plot!(af',linewidth=5.,color=colorant"lightgrey")
    plot!(x,v,marker=:x)
    plot!(x,D*v,marker=:o)
end
Out[9]:
In [14]:
a = -0.01;
y0 = 1.
N = 5
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*eye(N)) \(-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)
Out[14]:
In [ ]:
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(eye(2),D1)-kron(A,eye(N)))\(-kron(y0,d0))
U[1:end-1,:] = u2
U[end,:] = y0
plot(x,U)
In [ ]:
@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(eye(2),D1)-kron(A,eye(N)))\(-kron(y0,d0))
    U[1:end-1,1:end] = u2
    U[end,:] = y0
    plot(x,U)
    plot!(title=@sprintf("N=%i",N))
end
In [ ]:
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(eye(2),D1)-kron(A,eye(N)))\(-kron(y0,d0))
U[1:end-1,:] = u2
U[end,:] = y0
plot(x,U)
#[U[1,:]; (expm((A/25)*50)*y0)' ]
In [ ]:
N = 300
x = chebpts(N)
Ue = zeros(N+1,2)
A = [-0.01 1; -1 -0.01];
for i=1:N
    Ue[i,:] = (expm(A*(x[i]+1)*25)*y0)'
end
plot(x,sum(Ue.^2,2))
A