using Plots
using LinearAlgebra
function mgs(A)
m,n = size(A)
Q = copy(A)
R = zeros(n,n)
for i=1:n
R[i,i]=sqrt(transpose(Q[:,i])*Q[:,i])
if (R[i,i]<eps())
break
end
Q[:,i]=Q[:,i]/R[i,i]
for j=i+1:n
R[i,j]=transpose(Q[:,i])*A[:,j]
Q[:,j]=Q[:,j] - R[i,j] * Q[:,i]
end
end
return Q,R
end
mgs (generic function with 1 method)
function runge_f(t)
f_t = 1 ./ (25*(t .^2) .+1)
return f_t
end
runge_f (generic function with 1 method)
function monomial(n,m)
h = 2.0/(m-1)
t = [(i*h-1) for i=0:(m-1)]
A = ones(m,1)
for i=1:n-1
A = [A (t .^ i)]
end
return A
end
monomial (generic function with 1 method)
function monomial_coef(n,m)
h = 2.0/(m-1)
t = [(i*h-1) for i=0:(m-1)]
A = monomial(n,m)
y = runge_f.(t)
x = A\y
return x
end
monomial_coef (generic function with 1 method)
M = 512;
H = 2.0/(M-1);
T = [(i*H-1) for i=0:(M-1)];
Y = runge_f(T);
m_values = 2 .^ collect(4:8);
m_size = length(m_values);
err = Array{Float64}(undef, m_size);
for i=1:4
n = 2^(i+1)
m=96
h = 2.0/(m-1)
t = [(i*h-1) for i=0:(m-1)]
A = ones(m,1)
for i=1:n-1
A = [A (t .^ i)]
end
display(plot(A, legend=false))
end
n = 4;
B = monomial(n,M)
Approx = zeros(M,m_size)
for i=1:m_size
m = m_values[i]
x = monomial_coef(n,m)
Approx[:,i] = B*x
end
plot(T,Approx)
plot!(T,Y,label="Expected")
for i=1:m_size
err_abs = abs.(Y-Approx[:,i])
err[i] = dot(err_abs,err_abs)
end
err = log.(err .^0.5)
display(plot(log.(m_values),err, label="error"))
num = err[3:m_size] .- err[2:(m_size-1)];
den=err[2:m_size-1]-err[1:m_size-2];
q = sum(num ./ den)/(m_size-2);
s = exp(sum(err[2:m_size]-q*err[1:m_size-1])/(m_size-1));
print("The order of convergence is estimated to be")
print("\n")
print(q)
print("\n","\n")
print("The rate of convergence is")
print("\n")
print(s)
The order of convergence is estimated to be 0.2460158746597262 The rate of convergence is 3.035619860388405
n = 8;
B = monomial(n,M)
Approx = zeros(M,m_size)
for i=1:m_size
m = m_values[i]
x = monomial_coef(n,m)
Approx[:,i] = B*x
end
plot(T,Approx)
plot!(T,Y,label="Expected")
for i=1:m_size
err_abs = abs.(Y-Approx[:,i])
err[i] = dot(err_abs,err_abs)
end
err = log.(err .^0.5)
display(plot(log.(m_values),err, label="error"))
num = err[3:m_size] .- err[2:(m_size-1)];
den=err[2:m_size-1]-err[1:m_size-2];
q = sum(num ./ den)/(m_size-2);
s = exp(sum(err[2:m_size]-q*err[1:m_size-1])/(m_size-1));
print("The order of convergence is estimated to be")
print("\n")
print(q)
print("\n","\n")
print("The rate of convergence is")
print("\n")
print(s)
The order of convergence is estimated to be 0.40277320620291507 The rate of convergence is 1.5031719274633308
n = 16;
B = monomial(n,M)
Approx = zeros(M,m_size)
for i=1:m_size
m = m_values[i]
x = monomial_coef(n,m)
Approx[:,i] = B*x
end
plot(T,Approx)
plot!(T,Y,label="Expected")
for i=1:m_size
err_abs = abs.(Y-Approx[:,i])
err[i] = dot(err_abs,err_abs)
end
err = log.(err .^0.5)
display(plot(log.(m_values),err, label="error"))
num = err[3:m_size] .- err[2:(m_size-1)];
den=err[2:m_size-1]-err[1:m_size-2];
q = sum(num ./ den)/(m_size-2);
s = exp(sum(err[2:m_size]-q*err[1:m_size-1])/(m_size-1));
print("The order of convergence is estimated to be")
print("\n")
print(q)
print("\n","\n")
print("The rate of convergence is")
print("\n")
print(s)
The order of convergence is estimated to be -1.2579998394391712 The rate of convergence is 5.536127644117644
n = 32;
B = monomial(n,M)
Approx = zeros(M,m_size)
for i=1:m_size
m = m_values[i]
x = monomial_coef(n,m)
Approx[:,i] = B*x
end
plot(T,Approx)
plot!(T,Y,label="Expected")
for i=1:m_size
err_abs = abs.(Y-Approx[:,i])
err[i] = dot(err_abs,err_abs)
end
err = log.(err .^0.5)
display(plot(log.(m_values),err, label="error"))
num = err[3:m_size] .- err[2:(m_size-1)];
den=err[2:m_size-1]-err[1:m_size-2];
q = sum(num ./ den)/(m_size-2);
s = exp(sum(err[2:m_size]-q*err[1:m_size-1])/(m_size-1));
print("The order of convergence is estimated to be")
print("\n")
print(q)
print("\n","\n")
print("The rate of convergence is")
print("\n")
print(s)
The order of convergence is estimated to be -1.2579998394391712 The rate of convergence is 5.536127644117644
function bspline(n,m,M)
h_m = 2.0/(m-1)
h_M = 2.0/(M-1)
t = [(i*h_m-1) for i=0:(m+1)]
B = zeros(M,n)
for i=1:n
j_min = Int(max(ceil( (i-1)*h_m/h_M ),1))
j_med = Int(ceil( i*h_m/h_M )-1)
j_max = Int(min(ceil( (i+1)*h_m/h_M )-1,M))
for j=j_min:Int(min(j_med,M))
B[j,i] = (j*h_M-1-t[i])/h_m
end
if j_med+1 < M
for j=(j_med+1):j_max
B[j,i] = (t[i+2]-j*h_M+1)/h_m
end
end
end
return B
end
bspline (generic function with 1 method)
function bspline(n,m,M)
h_m = 2.0/(m-1)
h_M = 2.0/(M-1)
t = [(i*h_m-1) for i=0:(m+1)]
B = zeros(M,n)
for i=1:n
j_min = Int(max(ceil( (i-1)*h_m/h_M ),1))
j_med = Int(ceil( i*h_m/h_M )-1)
j_max = Int(min(ceil( (i+1)*h_m/h_M )-1,M))
for j=j_min:Int(min(j_med,M))
B[j,i] = (j*h_M-1-t[i])/h_m
end
if j_med+1 < M
for j=(j_med+1):j_max
B[j,i] = (t[i+2]-j*h_M+1)/h_m
end
end
end
return B
end
bspline (generic function with 1 method)
# an alternative way to define the bspline basis
function bspline_alt(n,m,M)
p = max(n,m)
h_m = 2.0/(m-1)
h_M = 2.0/(M-1)
t = [(i*h_m-1) for i=0:(p+1)]
T = [(i*h_M-1) for i=0:(M+1)]
B = zeros(M,n)
for k=1:n
for i=1:M
for j=1:m
if T[i+1] >= t[k] && T[i+1] < t[k+1]
B[i,k]=(T[i+1]-t[k])/h_m
elseif T[i+1] >= t[k+1] && T[i+1] < t[k+2]
B[i,k]=(t[k+2]-T[i+1])/h_m
end
end
end
end
return B
end
bspline_alt (generic function with 1 method)
function bspline_coef(n,m)
h = 2.0/(m-1)
t = [(i*h-1) for i=0:(m-1)]
A = bspline_alt(n,m,m)
y = runge_f.(t)
x = A\y
return x
end
bspline_coef (generic function with 1 method)
for i=1:4
n = 2^(i+1)
m=96
h = 2.0/(m-1)
t = [(i*h-1) for i=0:(m-1)]
A = bspline_alt(n,m,m)
display(plot(A, legend=false))
end
# testing our first function to define the basis
B = bspline(16,16,512)
plot(B)
# testing the alternative to define the absis
B = bspline_alt(16,16,512)
plot(B)
M = 512;
H = 2.0/(M-1);
T = [(i*H-1) for i=0:(M-1)];
Y = runge_f(T);
m_values = 2 .^ collect(4:8);
m_size = length(m_values);
err = Array{Float64}(undef, m_size);
n = 4;
Approx = zeros(M,m_size)
for i=1:m_size
m = m_values[i]
x = bspline_coef(n,m)
B = bspline_alt(n,m,M)
Approx[:,i] = B*x
end
plot(T,Approx)
plot!(T,Y,label="Expected")
for i=1:m_size
err_abs = abs.(Y-Approx[:,i])
err[i] = dot(err_abs,err_abs)
end
err = log.(err .^0.5)
display(plot(log.(m_values),err, label="error"))
num = err[3:m_size] .- err[2:(m_size-1)];
den=err[2:m_size-1]-err[1:m_size-2];
q = sum(num ./ den)/(m_size-2);
s = exp(sum(err[2:m_size]-q*err[1:m_size-1])/(m_size-1));
print("The order of convergence is estimated to be")
print("\n")
print(q)
print("\n","\n")
print("The rate of convergence is")
print("\n")
print(s)
The order of convergence is estimated to be 0.26876335831183384 The rate of convergence is 4.9648641695142
n = 8;
Approx = zeros(M,m_size)
for i=1:m_size
m = m_values[i]
x = bspline_coef(n,m)
B = bspline_alt(n,m,M)
Approx[:,i] = B*x
end
plot(T,Approx)
plot!(T,Y,label="Expected")
for i=1:m_size
err_abs = abs.(Y-Approx[:,i])
err[i] = dot(err_abs,err_abs)
end
err = log.(err .^0.5)
display(plot(log.(m_values),err, label="error"))
num = err[3:m_size] .- err[2:(m_size-1)];
den=err[2:m_size-1]-err[1:m_size-2];
q = sum(num ./ den)/(m_size-2);
s = exp(sum(err[2:m_size]-q*err[1:m_size-1])/(m_size-1));
print("The order of convergence is estimated to be")
print("\n")
print(q)
print("\n","\n")
print("The rate of convergence is")
print("\n")
print(s)
The order of convergence is estimated to be 0.15230817348724346 The rate of convergence is 6.583539724366618
n = 16;
Approx = zeros(M,m_size)
for i=1:m_size
m = m_values[i]
x = bspline_coef(n,m)
B = bspline_alt(n,m,M)
Approx[:,i] = B*x
end
plot(T,Approx)
plot!(T,Y,label="Expected")
for i=1:m_size
err_abs = abs.(Y-Approx[:,i])
err[i] = dot(err_abs,err_abs)
end
err = log.(err .^0.5)
display(plot(log.(m_values),err, label="error"))
num = err[3:m_size] .- err[2:(m_size-1)];
den=err[2:m_size-1]-err[1:m_size-2];
q = sum(num ./ den)/(m_size-2);
s = exp(sum(err[2:m_size]-q*err[1:m_size-1])/(m_size-1));
print("The order of convergence is estimated to be")
print("\n")
print(q)
print("\n","\n")
print("The rate of convergence is")
print("\n")
print(s)
The order of convergence is estimated to be 0.7302990696097617 The rate of convergence is 2.0342532903922033
n = 32;
Approx = zeros(M,m_size)
for i=1:m_size
m = m_values[i]
x = bspline_coef(n,m)
B = bspline_alt(n,m,M)
Approx[:,i] = B*x
end
plot(T,Approx)
plot!(T,Y,label="Expected")
for i=1:m_size
err_abs = abs.(Y-Approx[:,i])
err[i] = dot(err_abs,err_abs)
end
err = log.(err .^0.5)
display(plot(log.(m_values),err, label="error"))
num = err[3:m_size] .- err[2:(m_size-1)];
den=err[2:m_size-1]-err[1:m_size-2];
q = sum(num ./ den)/(m_size-2);
s = exp(sum(err[2:m_size]-q*err[1:m_size-1])/(m_size-1));
print("The order of convergence is estimated to be")
print("\n")
print(q)
print("\n","\n")
print("The rate of convergence is")
print("\n")
print(s)
The order of convergence is estimated to be -0.39470458104387757 The rate of convergence is 9.623611079944993
function legendre(n,m)
h = 2.0/(m-1)
t = [(i*h-1) for i=0:m-1]
A = monomial(n,m)
Q,R = qr(A)
S = diagm(1.0 ./ Q[m,:])
P = Q*S
return P[:,1:n]
end
legendre (generic function with 1 method)
function legendre_coef(n,m)
h = 2.0/(m-1)
t = [(i*h-1) for i=0:(m-1)]
A = legendre(n,m)
y = runge_f.(t)
x = A\y
return x
end
legendre_coef (generic function with 1 method)
M = 512;
H = 2.0/(M-1);
T = [(i*H-1) for i=0:(M-1)];
Y = runge_f(T);
m_values = 2 .^ collect(4:8);
m_size = length(m_values);
err = Array{Float64}(undef, m_size);
for i=1:4
n = 2^(i+1)
m=96
h = 2.0/(m-1)
t = [(i*h-1) for i=0:(m-1)]
P = legendre(n,m)
display(plot(P, legend=false))
end
n = 4;
B = legendre(n,M)
Approx = zeros(M,m_size)
for i=1:m_size
m = m_values[i]
x = legendre_coef(n,m)
Approx[:,i] = B*x
end
plot(T,Approx)
plot!(T,Y,label="Expected")
for i=1:m_size
err_abs = abs.(Y-Approx[:,i])
err[i] = dot(err_abs,err_abs)
end
err = log.(err .^0.5)
display(plot(log.(m_values),err, label="error"))
num = err[3:m_size] .- err[2:(m_size-1)];
den=err[2:m_size-1]-err[1:m_size-2];
q = sum(num ./ den)/(m_size-2);
s = exp(sum(err[2:m_size]-q*err[1:m_size-1])/(m_size-1));
print("The order of convergence is estimated to be")
print("\n")
print(q)
print("\n","\n")
print("The rate of convergence is")
print("\n")
print(s)
The order of convergence is estimated to be 0.2408324715806308 The rate of convergence is 3.058988236991122
n = 8;
B = legendre(n,M)
Approx = zeros(M,m_size)
for i=1:m_size
m = m_values[i]
x = legendre_coef(n,m)
Approx[:,i] = B*x
end
plot(T,Approx)
plot!(T,Y,label="Expected")
for i=1:m_size
err_abs = abs.(Y-Approx[:,i])
err[i] = dot(err_abs,err_abs)
end
err = log.(err .^0.5)
display(plot(log.(m_values),err, label="error"))
num = err[3:m_size] .- err[2:(m_size-1)];
den=err[2:m_size-1]-err[1:m_size-2];
q = sum(num ./ den)/(m_size-2);
s = exp(sum(err[2:m_size]-q*err[1:m_size-1])/(m_size-1));
print("The order of convergence is estimated to be")
print("\n")
print(q)
print("\n","\n")
print("The rate of convergence is")
print("\n")
print(s)
The order of convergence is estimated to be 0.38479479155555435 The rate of convergence is 1.521563387971189
n = 16;
B = legendre(n,M)
Approx = zeros(M,m_size)
for i=1:m_size
m = m_values[i]
x = legendre_coef(n,m)
Approx[:,i] = B*x
end
plot(T,Approx)
plot!(T,Y,label="Expected")
for i=1:m_size
err_abs = abs.(Y-Approx[:,i])
err[i] = dot(err_abs,err_abs)
end
err = log.(err .^0.5)
display(plot(log.(m_values),err, label="error"))
num = err[3:m_size] .- err[2:(m_size-1)];
den=err[2:m_size-1]-err[1:m_size-2];
q = sum(num ./ den)/(m_size-2);
s = exp(sum(err[2:m_size]-q*err[1:m_size-1])/(m_size-1));
print("The order of convergence is estimated to be")
print("\n")
print(q)
print("\n","\n")
print("The rate of convergence is")
print("\n")
print(s)
The order of convergence is estimated to be 1.0650633160226346 The rate of convergence is 0.6086171218702784
n = 32;
B = legendre(n,M)
Approx = zeros(M,m_size-1)
for i=1:m_size-1
m = m_values[i+1]
x = legendre_coef(n,m)
Approx[:,i] = B*x
end
plot(T,Approx)
plot!(T,Y,label="Expected")
err = zeros(m_size-1)
for i=1:m_size-1
err_abs = abs.(Y-Approx[:,i])
err[i] = dot(err_abs,err_abs)
end
err = log.(err .^ 0.5)
a = [4; 4.5; 5]
plot(log.(m_values[1:m_size-1]),err, label="error")
num = err[3:m_size-1] .- err[2:(m_size-2)];
den=err[2:m_size-2]-err[1:m_size-3];
q = sum(num ./ den)/(m_size-3);
s = exp(sum(err[2:m_size-1]-q*err[1:m_size-2])/(m_size-2));
print("The order of convergence is estimated to be")
print("\n")
print(q)
print("\n","\n")
print("The rate of convergence is")
print("\n")
print(s)
display(plot!(a, (a*(-1.5) .+ 6), label="slope=-1.5"))
The order of convergence is estimated to be 1.396601993488745 The rate of convergence is 0.46366810375093176