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