assignment 5.ipynb

Kernel: Python 3 (Anaconda)

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt

n = 100
m = 10

# Define G
G = 0.5*np.diag(np.ones(n)) + 0.25*np.diag(np.ones(n-1),k=-1) + 0.25*np.diag(np.ones(n-1),k=1)

# Define K
K = np.linalg.matrix_power(G,m)

# define Heaviside
H = lambda x: np.piecewise(x, [x < 0, x >= 0], [0, 1])

# define u
t = np.arange(1,n+1)/(n+1)
u1 = H(t-0.4) - H(t-0.6)
u2 = np.exp(-100*(t-0.5)**2)

# noisy data
sigma = 0
noise = sigma*np.random.randn(n)
f1 = K@u1
f2 = K@u2

# plot
plt.subplot(121)
plt.plot(t,u1,t,f1)

plt.subplot(122)
plt.plot(t,u2,t,f2)

[<matplotlib.lines.Line2D at 0x7fa63043a940>,
 <matplotlib.lines.Line2D at 0x7fa6303dd9e8>]

U, S, V = np.linalg.svd(K)
# noise
plt.semilogy(S,'*')
plt.semilogy(U.T@f2,'o')

[<matplotlib.lines.Line2D at 0x7fa630075198>]

alpha = 1e-16
F = (V[0:n,:][e-mail](S/(S**2 + alpha))@U[:,0:n].T)
plt.plot(t,u1,t,F@f1)

[<matplotlib.lines.Line2D at 0x7fa63005dc18>,
 <matplotlib.lines.Line2D at 0x7fa63005dda0>]

alpha = 1e-16
F = (V[0:n,:][e-mail](S/(S**2 + alpha))@U[:,0:n].T)
plt.plot(t,u2,t,F@f2)

[<matplotlib.lines.Line2D at 0x7fa6301ae780>,
 <matplotlib.lines.Line2D at 0x7fa6301ae908>]