Matlab code

Graph of p-norm equals to one

$\Vert x{\Vert }_p=1,x\in {\mathbb{R}}^2$

%% The graph of p-norm equals to one
%
p=4;                      % define the value of p
theta = linspace(0, pi/2);
x = cos(theta).^(2/p);
y = sin(theta).^(2/p);
plot(x, y, -x, y, x, -y, -x, -y, 'LineWidth', 2)
axis image

function approximation

Fourier series of

$f(x)=x$

%%
% Fourier series approximation of f(x)=x on -pi<x<pi
%

N = 10^4;                 % number of sampling points
x = linspace(-pi, pi, N); % sampling points
f = x;                    % f(x)=x
g = zeros(1, N);
M = 10;                   % Approximate by M term Fourier series
error_2 = zeros(1, M);
error_inf = zeros(1, M);
for kk=1:M
    g = g-2*sin(kk*x)*(-1)^kk/kk;
    error_2(kk) = norm(g-f)*2*pi/N;
    error_inf(kk) = norm(g-f, 'Inf');
end

%% compare the function and its approximation
plot(x, f, x, g, 'linewidth', 1)
title([num2str(M) ' term Fourier series'])
legend('f', 'Fourier series')
xlabel('x')
axis image

%% plot the errors
figure; 
plot(1:M, error_2, 1:M, error_inf, 'linewidth', 2)
legend('error in 2-norm', 'error in Inf-norm')
xlabel('number of terms')

Fourier series of

$f(x)=x^2$

• You can change
  $f(x)$ to see the Fourier approximation of arbitrary function

%%
% Fourier series approximation of f(x)=x^2 on -pi<x<pi using FFT
%

N = 10^4;                 % number of sampling points
x = linspace(-pi, pi, N+1); x = x(1:N); % sampling points
f = x.^2;                 % f = x^2
M = 3;                   % Approximate by M term Fourier series

%% Fourier approximation
fft_f = fft(f);
fft_f(M+1:N-M+1) = 0;
g = real(ifft(fft_f));

%% compare the function and its approximation
plot(x, f, x, g, 'linewidth', 1)
title([num2str(M) ' term Fourier series'])
legend('f', 'Fourier series')
xlabel('x')

Timing for

$Ax=b$

In the following we compare the time difference between perform

$A\text{b}$ and

$A^{-1}b$

clear;
clc;
n = 3000;                    % size of the linear system
A = rand(n); b = rand(n,1);
m=10;                        % number of trial

tic
for ii=1:m
    x = A\b;
end
disp(['The average time for x=A\b is  ', num2str(toc/m)])

tic
for ii=1:m
    x = inv(A)*b;
end
disp(['The average time for x=inv(A)*b is  ', num2str(toc/m)])

Multidimensional scaling

Example in

${\mathbb{R}}^2$ with random sampling

%% 資料生成
% 在單位圓上任取十點做成一多邊形
n = 10;
theta = sort(rand(n,1))'*2*pi;
X = [cos(theta); sin(theta)];
X = X + rand(2,1)*ones(1,n);
plot(X(1,:), X(2,:), 'o-', 'linewidth', 1); 
axis image; 
hold on;

%% EDM
D = squareform(pdist(X', 'squaredeuclidean'));  % 求出其 EDM

%%
% n 個 samples
mu = sum(D, 1)/n; D = D - mu;           %
mu = sum(D, 2)/n; B = D - mu;           %
B = -0.5*B;                             % B = -0.5*H*D*H
[V, S] = eig(B, 'vector');              % 求出特徵值及特徵向量
[sqD, ind] = sort(sqrt(S), 'descend');  % 將特徵值按大小排列
sqD = sqD(1:2);                         % 取前 2 個
V = V(:, ind(1:2));                     % 取相對應的特徵向量
Y = V'.*(sqD*ones(1,n));                % 求出 k 維資料點

plot(Y(1,:), Y(2,:), 'x-', 'linewidth', 1);
legend('original data', 'MDS recovered data')

Linear transformation

%% 
% Define the matrix of linear transformation from R^2 to R^2

A = [1 -2; 3 2]

%% construct an unit circle
%
x = linspace(0, pi/2);
x1 = [cos(x); sin(x)];
x2 = [cos(x+pi/2); sin(x+pi/2)];
x3 = [cos(x+pi); sin(x+pi)];
x4 = [cos(x+pi*3/2); sin(x+pi*3/2)];
subplot(1,2,1)
plot(x1(1,:), x1(2,:), 'b', x2(1,:), x2(2,:), 'r', x3(1,:), x3(2,:), 'g', x4(1,:), x4(2,:), 'k', 'linewidth', 2)
hold on
axis equal

%% construct the vertical and horizontal lines
y = linspace(-1, 1, 6);
for ii=1:6
    plot(y(ii)*[1 1], [-1 1], 'b')
    plot([-1 1], y(ii)*[1 1], 'r')
end

%% construct the transformed unit circle
z1 = A*x1;
z2 = A*x2;
z3 = A*x3;
z4 = A*x4;
subplot(1,2,2)
plot(z1(1,:), z1(2,:), 'b', z2(1,:), z2(2,:), 'r', z3(1,:), z3(2,:), 'g', z4(1,:), z4(2,:), 'k', 'linewidth', 2)
hold on
axis equal

%% construct the transformed vertical and horizontal lines
for ii=1:6
    y1 = A*[y(ii)*[1 1]; [-1 1]];
    plot(y1(1,:), y1(2,:), 'b')
    y2 = A*[[-1 1]; y(ii)*[1 1]];
    plot(y2(1,:), y2(2,:), 'r')
end

%% Not yet finished
%
%% 
% Define the matrix of linear transformation from R^2 to R^2

A = [1 0; 0.1 1.1]

%% construct an unit circle
%
x = linspace(0, 2*pi);
x1 = [cos(x); sin(x)];
plot(x1(1,:), x1(2,:), 'b', 'linewidth', 1)
hold on

%% consider trajectory of four points
x = zeros(2, 101);
M = 100;
for jj=0:M
    x(:,1) = [cos(2*pi*jj/M); sin(2*pi*jj/M)];
    for ii=1:100
        x(:,ii+1) = A*x(:,ii);
    end
    plot(x(1,:), x(2,:), 'linewidth', 1)
end

axis([-10 10 -10 10])
axis equal

A = [1 0; 3 2]

%% construct an unit circle
%
x = linspace(0, pi/2);
x1 = [cos(x); sin(x)];
x2 = [cos(x+pi/2); sin(x+pi/2)];
x3 = [cos(x+pi); sin(x+pi)];
x4 = [cos(x+pi*3/2); sin(x+pi*3/2)];
subplot(1,2,1)
plot(x1(1,:), x1(2,:), 'b', x2(1,:), x2(2,:), 'r', x3(1,:), x3(2,:), 'g', x4(1,:), x4(2,:), 'k', 'linewidth', 2)
hold on
axis equal

[V,D] = eig(A, 'vector'); 

if isreal(D(1,1))==0
    %% construct the vertical and     horizontal lines
    y = linspace(-1, 1, 6);
    for ii=1:6
        plot(y(ii)*[1 1], [-1 1], 'b')
        plot([-1 1], y(ii)*[1 1], 'r')
    end
else
    %% construct the vertical and horizontal lines
    y1 = V(:,1)*[1 -1];
    y1 = y1/norm(y1);
    plot(y1(1,:), y1(2,:), 'b')
    y2 = V(:,2)*[1 -1];
    y2 = y2/norm(y2);
    plot(y2(1,:), y2(2,:), 'r')
end

%% construct the transformed unit circle
z1 = A*x1;
z2 = A*x2;
z3 = A*x3;
z4 = A*x4;
subplot(1,2,2)
plot(z1(1,:), z1(2,:), 'b', z2(1,:), z2(2,:), 'r', z3(1,:), z3(2,:), 'g', z4(1,:), z4(2,:), 'k', 'linewidth', 2)
hold on
axis equal

%% construct the transformed vertical and horizontal lines
if isreal(D(1,1))==0
    %% construct the vertical and     horizontal lines
    y = linspace(-1, 1, 6);
    for ii=1:6
        plot(y(ii)*[1 1], [-1 1], 'b')
        plot([-1 1], y(ii)*[1 1], 'r')
    end
else
    %% construct the vertical and horizontal lines
    y1 = V(:,1)*[1 -1];
    y1 = A*y1/norm(y1);
    plot(y1(1,:), y1(2,:), 'b')
    y2 = V(:,2)*[1 -1];
    y2 = A*y2/norm(y2);
    plot(y2(1,:), y2(2,:), 'r')
end