%% Matlab codes for the course Applied numerical methods with Matlab (2019/20 spring)

%% Problem 13

C=[	2 -1 1 1 1
	-1 2 -4 2 6
	-3 3  2 -5 4
	 6 7 -1 4 -9
	 4 -3 0 -3 3]

%% Problem 14

format long

A=rand(6);
matrixnorm=norm(A)
spectral_radius=max(abs(eig(A)))

%% Problem 15

A=[0.5, 0.5; -0.5, 0.25]

%% Problem 16

clear x, clear y
x(1)=1; for i=1:40 x(i+1)=2-2*(2-x(i))/(4-x(i)^2); end
disp('iteration steps of the sequence x_i')
disp(x')
log_x_error=log(abs(x-sqrt(2)));

y(1)=1; y(2)=2; for i=2:15 y(i+1)=y(i)-(y(i)^2-2)/(y(i)+y(i-1)); end
disp('iteration steps of the sequence y_i')
disp(y')
log_y_error=log(abs(y-sqrt(2)));

plot(log_x_error(1:end-1),log_x_error(2:end),log_y_error(1:end-1),log_y_error(2:end))
axis equal

legend('x','y')

%% Problem 17

d=3.4e-16;
a=1+d;
b=sqrt(a)
x=b-1

%% Problem 19

x=-25;
e=0;
for i=0:200
    e=e+x^i/factorial(i);
end
e

%% Problem 20

x=10^8;
v=1/(sqrt(x^2+1)-x)

%% Problem 21

x=10^-15;
v=((1+x)-1)/x

%% Problem 24

A=toeplitz([2,-1,0,0]);

%% Problem 28

A=toeplitz([3,-1,-1,0]);
b=[1;1;1;1];
Atilde=A+(0.04*rand(4)-0.02);
x=A\b
xtilde=Atilde\b
rel_error=norm(x-xtilde,inf)/norm(x,inf)

%% Problem 29

A=[34,55;55,89];
b=[21;34];
x1=[-0.11;0.45];
x2=[-0.99,1.01];

%% Problem 30

A=[1,1,1;2,4,4;-1,5,8];
b=[4,14,17];

% Code for the LU decomposition
% The code rewrites A with L and U

[n,n] = size(A);
tol = 1.e-6;

for k = 1:n-1
   if abs(A(k,k)) < tol
      disp(['Small pivot encountered in column ' int2str(k)])
   end
   for i = k+1:n
      A(i,k) = A(i,k)/A(k,k);
      A(i,k+1:n) = A(i,k+1:n) - A(i,k)*A(k,k+1:n);
   end
end

U=triu(A)
L=eye(n)+tril(A,-1)

%% Problem 32

format rat
A=[1/10000,1;1,1]
[L,U,P]=lu(A)

%% Problem 33

format rat
A=[4,2,1;2,17,13/2;1,13/2,23/2];
lu(A)
chol(A)

%% Problem 35

format rat
A=[4,2,1;2,17,13/2;1,13/2,23/2];
b=[3;-9/2;-21/2];
G=chol(A); G=G'
y=G\b
x=G'\y

%% Problem 36

format long
A=[-1,5,-2;1,1,-4;4,-1,2]
b=[5;1;-1]

L=-tril(A,-1)
R=-triu(A,1)
D=diag(diag(A))

x=zeros(3,1); % initial vector
for i=1:25
    pause(0.3)
    i
    x=x % modify according to the method
end

%% Problem 37

% Computation of the optimal omega parameter

A=[2,-1;-1,2], b=[1;3]

syms omega

L=-tril(A,-1)
R=-triu(A,1)
D=diag(diag(A))

B=(1/omega*D)\(L+R+(1-omega)/omega*D)
f=((1/omega)*D)\b
eigenvalues=eig(B) % Eigenvalues of the iteration matrix

o=0:0.05:2;
plot(o,max(abs([1-o/2;1-3*o/2]))) % graph of the spectral radius against omega

%%

% Checking the convergence speed

x_exact=A\b;

omega=1
Bomega=eval(B);
fomega=eval(f);

x=zeros(2,1); % initial vector

i=0;
while norm(x-x_exact)>=10^-6
    i=i+1
    x=Bomega*x+fomega 
end

eigenvalues_omega=eval(eigenvalues)

%% Problem 38

x=[2;6;-3];
v=x; v(1)=v(1)+norm(v);
H=eye(3)-2*v*v'/(v'*v)
H*x

c=6/norm(x(2:3)); s=3/norm(x(2:3));
G=blkdiag(1,[c,-s;s,c])
G*x

%% Problem 39

A=[0,0;1,3;0,2]
x=A(1:2,1);
c=x(1)/norm(x), s=-x(2)/norm(x)
G1=blkdiag([c,-s;s,c],1) % The first Givens matrix
A1=G1*A
x=A1(2:3,2);
c=x(1)/norm(x), s=-x(2)/norm(x)
G2=blkdiag(1,[c,-s;s,c]) % The second Givens matrix
R=G2*A1 % R
Q=G1'*G2' % Q
Q*R % Check whether A=QR 
qr(A)

%% Problem 41

A=[0,0;1,3;0,2] % the matrix from the previous problem
b=[3;4;1];

% mldivide

A\b

% normal equation

(A'*A)\(A'*b)

% QR

[Q,R]=qr(A);
c=Q'*b;
R(1:2,1:2)\c(1:2)

%% Problem 42

A=[1,1,1;
   0,0,1;
   1,1,1;
   4,2,1]

b=[1;0;3;2]

A\b %-x^2+3x

% with polynomial fitting

z=polyfit([1,0,1,2],[1,0,3,2],2)

%% Problem 43

A=[4.2,-0.1,0.2,0.1;0.05,3,-0.1,-0.05;0.5,-0.5,2,0.1;0.1,0.2,0.3,-1]
eig(A)

%% Problem 44

% a)

A=[3,1;1,3]
eig(A) % eigenvalues of A
A(2,1)=A(2,1)+0.1
eig(A) % eigenvalues of the perturbed matrix

% Change of eigenvalues
change_of_eigenvalues=norm([0,0;0.1,0])


%%

% b) 

A=[9,2;1,3]
eig(A) % eigenvalues of A

A(2,1)=A(2,1)+0.1
eig(A) % eigenvalues of the perturbed matrix

[V,S]=eig(A)
inv(V)*A*V % A can be diagonalized with the matrix of the eigenvectors.

% Change of eigenvalues
change_of_eigenvalues=norm([0,0;0.1,0])*cond(V)

%% Problem 45

A = [2,1,0;1,2,1;0,1,2]
eig(A)

%%
% power method

y=rand(size(A,1),1);

for i=1:50
  i
  x=A*y;
  y=x/norm(x)
  nu=y'*A*y
  pause
end

%% Problem 46

A=toeplitz([20:-1:1]);

ev=eig(A)

tic
for i=1:100
    [Q,R]=qr(A);
    A=R*Q
    pause
end
toc

ev=diag(A)
ev_error=sum(abs(A-diag(diag(A))))'

%% Problem 47

A=[0,1;1,0];

for i=1:100
    [Q,R]=qr(A);
    A=R*Q
    pause
end

%% Problem 49

format long
xstar=((27^(1/2)*5240^(1/2))/27 + 352/27)^(1/3) - 26/(9*((27^(1/2)*5240^(1/2))/27 + 352/27)^(1/3)) - 2/3

syms x
p=x^3+2*x^2+10*x-20
pder=diff(p)
p=matlabFunction(p);
pder=matlabFunction(pder);
fplot(p,[1,2])

%%
% b)
x=2, 
for i=1:6 
    i
    x=x-p(x)/pder(x) 
    disp('Press <enter>')
    pause
end

%%
% c)
x=2, for i=1:50 i, x=20/(x^2+2*x+10), disp('Press <enter>'), pause, end

%%
% d)

fzero(p,1,optimset('Display','iter'))
pause

fzero(p,[1,2],optimset('Display','iter'))

%% Problem 50

clear error
syms x
f=x^3-5*x^2+8*x-4;
fder=diff(f);
f=matlabFunction(f);
fder=matlabFunction(fder);

x=3, for i=1:15 i, x=x-f(x)/fder(x), error(i)=abs(x-2); disp('Press <enter>'), pause, end

plot(log(error(1:end-1)),log(error(2:end)))
axis equal

%% Problem 51

f=@(x)x.*(63*x.^4-70*x.^2+15)/8;
F=@(x)x-f(x)/10;

Fs=sym(F);
Fders=diff(Fs); Fder=matlabFunction(Fders);

fplot(F,[0.8,1]), axis equal, title('Graph of F'), hold on, fplot(@(x)x,[0.8,1])
pause, hold off
fplot(Fder,[0.8,1]), title('Graph of F''')
q=Fder(0.8)

x0=1, x1=F(x0)
k_est=log(5e-7/abs(x1-x0)*(1-q))/log(q)

pause

x=1, for k=1:112 k, x=F(x), end

%% Problem 52

xstar=2;
F=@(x)x-1+1/(x-1);
Fs=sym(F)
Fders=diff(Fs)
Fder=matlabFunction(Fders);
Fder(xstar) % This shows at least second order convergence

pause

fplot(F,[1.5,2.5]), hold on, fplot(@(x)x,[1.5,2.5]), title('Graph of F')
hold off

%% Problem 53

syms x y z 
Fs=[x^2+2*y^2-y-2*z;x^2-8*y^2+10*z;x^2-7*y*z]
Js=jacobian(Fs)
F=matlabFunction(Fs,'vars',{[x;y;z]})
J=matlabFunction(Js,'vars',{[x;y;z]})

x=[1;1;1],
for k=1:8 z=J(x)\F(x); x=x-z, end

pause

fsolve(F,[1;1;1],optimset('Display','iter'))

%% Problem 54

F=@(x)(1-x(1))^2+5*(x(2)-x(1)^2)^2
fminunc(F,[2,2])

%% Problem 55

x=[-1,0,1];
y=[6,3,2];
z=polyfit(x,y,2)
xx=linspace(min(x),max(x),500);
yy=polyval(z,xx);
plot(x,y,'ro',xx,yy,'b-')

%% Problem 56

f=@(x)log(x+1);
x0=0.45;
x=[0,0.6,0.9];
xx=linspace(0,0.9,100);
z=polyfit(x,f(x),2);
yy=polyval(z,xx);
plot(x,f(x),'or',xx,yy,'b-',xx,f(xx),'-k')
legend('nodes','interpolation polynomial','f(x)')
error_at_x0=(f(x0)-polyval(z,x0))

%% Problem 57

f=@(x)sin(x);
x0=0.3; % a given point where we calculate the interpolation error
n=4 % number of the sub-intervals
x=linspace(0,pi,n+1);
xx=linspace(0,pi,100);
yy=interp1(x,f(x),xx,'linear');
plot(x,f(x),'or',xx,yy,'b-',xx,f(xx),'-k')
legend('given points','interpolation polynomial','f(x)')
error_at_x0=f(x0)-interp1(x,f(x),x0,'linear')

%% Problem 58

np1=12; 
x0=0.3; % an arbitrary point where the error is checked
f=@(x)1./(1+x.^2);
x=cos(pi/2/np1+[0:np1-1]*pi/np1);
y=f(x);
z=polyfit(x,y,np1-1);
xx=linspace(-1,1,200);
yy=polyval(z,xx);
error_at_x0=abs(f(x0)-polyval(z,x0))

plot(x,y,'ro',xx,yy,'b-',xx,f(xx),'k-')
legend('given points','interpolation polynomial','f(x)')
pause

fs=sym(f)
fders=diff(fs,np1)
fder=matlabFunction(fders)
xx=linspace(-1,1,500);
plot(xx,fder(xx))

%% Problem 59

x=[1,3]; y=[2,1];
s=spline(x,[1,y,2]) % the first and the last elements in the second vector are the derivatives at the endpoints
s.coefs % the coefficients of the polynomial in (x-1) - local coordinates

xx=linspace(1,3,100);
plot(xx,ppval(s,xx),'b-',x,y,'ro')
axis equal

%% Problem 60

x=[1,2,3]; y=[2,1,1];
d=(1/3*[2,1,0;1,4,1;0,1,2])\[-1;-1;0]
s=spline(x,[d(1),y,d(end)]) % 
s.coefs 

xx=linspace(1,3,100);
plot(xx,ppval(s,xx),'b-',x,y,'ro')
axis equal

%% Problem 61

n=3;
x=[0:n]*2*pi/(n+1);
f=abs(x-pi);

plot(x,f,'o')

a0=sym(sum(f)/(n+1))
a1=sym(f*cos(x)'*2/(n+1))
a2=sym(f*cos(2*x)'*2/(n+1))

b1=sym(f*sin(x)'*2/(n+1))

t=@(x)a0+a1*cos(x)+b1*sin(x)+a2*cos(2*x);

hold on
fplot(t,[0,2*pi])
hold off

%% Problem 62

clear all
n=5;
f=[0,1,1,0,-1,-1];
x=sym([0:n]*2*pi/(n+1));
plot(x,f,'o')

%%%%%%%%%%%%%%%%%%%
% a)
% The real discrete Fourier coefficients can be calculated as follows:

a0=sym(sum(f)/(n+1))
a1=sym(f*cos(x)'*2/(n+1))
a2=sym(f*cos(2*x)'*2/(n+1))
a3=sym(f*cos(3*x)'/(n+1))

b1=sym(f*sin(x)'*2/(n+1))  % We need only this.
b2=sym(f*sin(2*x)'*2/(n+1))

t=@(x)a0+a1*cos(x)+b1*sin(x)+a2*cos(2*x)+b2*sin(2*x)+a3*cos(3*x);

hold on
fplot(t,[0,2*pi])
hold off

%%%%%%%%%%%%%%%%%%%
% b) 

c=sym(fft(f')/(n+1))

a0=sym(c(1))
a1=sym(2*real(c(2)))
a2=sym(2*real(c(3)))
a3=sym(c(4))

b1=sym(-2*imag(c(2)))
b2=sym(-2*imag(c(3)))

%%%%%%%%%%%%%%%%%%%
% c) 

w=exp(-1i*2*pi/(n+1));
[xi,yi]=meshgrid(0:n,0:n);
F=w.^(xi.*yi)
c=F*f'/(n+1) % computation of c with multiplication with the Fourier matrix

F_small=F([1,2,3],[1,3,5])
D=diag(w.^[0:2])
f1=f([1,3,5])'
f2=f([2,4,6])'

F_small_f1=F_small*f1;
DF_small_f2=D*(F_small*f2);

c1=(F_small_f1+DF_small_f2)/(n+1)
c2=(F_small_f1-DF_small_f2)/(n+1)


%% Problem 63

x=linspace(0,2*pi,257);
x=x(1:256);
f=sin(3*x).*cos(4*x)-sin(x).^2;
f=f+rand(1,256)*0.1;
f=f+rand(1,256)*0.3; % sample from a noisy signal
plot(x,f)
xlabel('time')
ylabel('signal value')
pause

f_hat=fft(f);
plot((0:19)/2/pi,abs(f_hat(1:20)))

% f_hat(12:246)=zeros(1,size([12:128],1));
% f_hat(130:246)=zeros(1,size([12:128],1));
% 
% ff=ifft(f_hat)
% 
% plot(x,ff)

%% Problem 64

x0=0.05;
h=[0.0016,0.0008];
f=@(x)1./x;

d=(f(x0+h)-f(x0+0*h))./h
error=abs(-400-d) % -400 is the exact derivative

rich=2*d(2)-d(1) 


%% Problem 65

x0=0.05;
h=[0.0016,0.0008];
f=@(x)1./x;

d=(-3*f(x0+0*h)+4*f(x0+h)-f(x0+2*h))./h/2
error=abs(-400-d) % -400 is the exact derivative

%% Problem 66

syms x
l_0=(x-1/4)*(x-1/2)*(x-3/4)*(x-1)/((0-1/4)*(0-1/2)*(0-3/4)*(0-1));
N40=int(l_0,0,1)
N44=N40
l_1=(x-0)*(x-1/2)*(x-3/4)*(x-1)/((1/4-0)*(1/4-1/2)*(1/4-3/4)*(1/4-1));
N41=int(l_1,0,1)
N43=N41
N42=(1-2*(N40+N41))

I=eval((pi-0)*(N40*sin(0)+N41*sin(pi/4)+N42*sin(pi/2)+N43*sin(3*pi/4)+N44*sin(pi)))

integral(@(x)sin(x),0,pi)

%% Problem 67

format long
n=2;
h=pi/n;
x=0:h/2:pi;
f=sin(x);

Iexact=2;
I=h/6*(f(1)+f(end)+4*sum(f(2:2:end-1))+2*sum(f(3:2:end-2)))
error=I-Iexact

format short

%% Problem 68

format long
x=cos(pi/8+[0:3]*pi/4);
syms t
f=(4*t^7+4*t^6-3*t^5+6*t^4-t^3+4*t^2-t+5); % an arbitrary 7th degree polynomial 
sf=f/sqrt(1-t^2);
I_exact=eval(int(sf,-1,1)) % the exact integral value 

fm=matlabFunction(f);
I=sum(fm(x))*pi/4
format short

%% Problem 69

syms x

Legnodes=solve(5*x^3-3*x) % These are the nodes

w=prod(x-Legnodes); % The nodal polynomial

for i=1:3 % computation of the weights (these are now given in the problem)
weight(i)=int(w/(x-Legnodes(i))/subs(diff(w),Legnodes(i)),-1,1);
end

weight

f=@(x)exp(-x.^2);
I_exact=integral(f,0,1) % The "exact" value must be computed with Matlab

n=8 % number of subintervals

h=1/n; % 1 is the lenght of the interval

I=0;
for i=1:3
I=I+sum(f(h/2*(1+eval(Legnodes(i))):h:1))*eval(weight(i));
end
I=I*h/2;

error=I_exact-I


%% Problem 70

f=@(x,y)atan(y);
x0=0, y0=1, xmax=1, N=20

h=(xmax-x0)/N;

x=zeros(1,N+1); y=x;
x(1)=x0;
y(1)=y0;

for i=1:N    % The EE iteration
    x(i+1)=x(i)+h;
    y(i+1)=y(i)+h*f(x(i),y(i));
end

plot(x,y,'bo-')

%% Problem 71

yexact=@(x)4*exp(x.^3/3)-1;
f=@(x,y)x.^2.*(1+y);
x0=0, y0=3,
xmax=1, N=4, xx=linspace(x0,xmax,500);

h=(xmax-x0)/N;

x=zeros(1,N+1); y=x;
x(1)=x0;
y(1)=y0;

for i=1:N  % The RK4 iteration
    x(i+1)=x(i)+h;
    k1=f(x(i),y(i)); k2=f(x(i)+h/2,y(i)+k1*h/2); k3=f(x(i)+h/2,y(i)+k2*h/2); k4=f(x(i)+h,y(i)+k3*h);
    y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);
    
    %The Heun iteration
    %k1=f(x(i),y(i)); k2=f(x(i)+h,y(i)+k1*h);       
    %y(i+1)=y(i)+h/2*(k1+k2);
end

plot(x,y,'bo-',xx,yexact(xx),'r-')

error=y(end)-yexact(xmax) % Error at xmax

%% Problem 72

f=@(x,y)[-0.013*y(1)-1000*y(1)*y(3);-2500*y(2)*y(3);-0.013*y(1)-1000*y(1)*y(3)-2500*y(2)*y(3)];

% RK4

clear x y
N=100000;
h=50/N;
x(1)=0;
y=[1;1;0];
for i=1:N
 x(i+1)=x(i)+h;
 k1=f(x(i),y(:,i)); k2=f(x(i)+h/2,y(:,i)+k1*h/2); k3=f(x(i)+h/2,y(:,i)+k2*h/2); k4=f(x(i)+h,y(:,i)+k3*h);
 y(:,i+1)=y(:,i)+h/6*(k1+2*k2+2*k3+k4);
end
plot(x,y)

pause

ode45(f,[0,50],[1,1,0])

pause

ode23s(f,[0,50],[1,1,0])

pause

% The implicit Euler method

N=50;
h=50/N;

x=zeros(1,N+1); 
x(1)=0;
y=zeros(3,N+1);
y(:,1)=[1;1;0];

J=@(y)h*[-0.013-1000*y(3)-1/h,0,-1000*y(1);0,-2500*y(3)-1/h,-2500*y(2);-0.013-1000*y(3),-2500*y(3),-1000*y(1)-2500*y(2)-1/h];

for i=1:N
    x(i+1)=x(i)+h;
    z=y(:,i);
    for k=1:10  % We solve the nonlinear system with the Newton method using 10 iteration steps
    z=z-J(z)\(y(:,i)+h*f(x(i+1),z)-z);
    end
    y(:,i+1)=z;
end

plot(x,y,'bo-')

%% Problem 73

% Checking the consistency order of a multistep method
%%%%%%%%%%%%%%%%%%

a=sym([3/2,-2,1/2]), b=sym([1,0,0]) % The coefficients of the method

steps=size(a,2)-1; % Number of steps
i=sym(steps:-1:0);

d=sum(a); k=1;
while d==0
    d=sum((i.^k).*a/factorial(k)-(i.^(k-1)).*b/factorial(k-1));
    k=k+1;
end

order=k-2

%%%%%%%%%%%%%%%%%

format long

yexact=@(x)4*exp(x.^3/3)-1;
f=@(x,y)x.^2.*(1+y);

x0=0, y0=3,
xmax=1, N=40, xx=linspace(x0,xmax,500);

h=(xmax-x0)/N;

x=zeros(1,N+1); y=x;
x(1)=x0;
y(1)=y0;
x(2)=x(1)+h;
y(2)=y(1)+h*f(x(1)+h/2,y(1)+h/2*f(x(1),y(1))); % RK2

for i=2:N
    x(i+1)=x(i)+h;
    z=y(i);
    error=1;
    while error>10^-6  % solution of the nonlinear equation   
    znew=2*h/3*f(x(i+1),z)+4/3*y(i)-1/3*y(i-1);
    error=abs(z-znew);
    z=znew;
    end
    y(i+1)=z;
end

final_error=y(end)-yexact(xmax)

plot(x,y,'bo-',xx,yexact(xx),'r-')

format short

%% Problem 74

yexact=@(x)4*exp(x.^3/3)-1;
f=@(x,y)x.^2.*(1+y);
x0=0, y0=3,
xmax=1, N=40, xx=linspace(x0,xmax,500);

h=(xmax-x0)/N;

x=zeros(1,N+1); y=x;
x(1)=x0;
y(1)=y0;
x(2)=x(1)+h;
y(2)=y(1)+h*f(x(1)+h/2,y(1)+h/2*f(x(1),y(1))); % RK2

for i=2:N
    x(i+1)=x(i)+h;
    y(i+1)=y(i)+h*(3/2*f(x(i),y(i))-1/2*f(x(i-1),y(i-1)));
end

final_error=y(end)-yexact(xmax)

plot(x,y,'bo-',xx,yexact(xx),'r-')

%% Problem 75

yexact=@(x)(2*x+1).*exp(x);
a=0; A=1; b=1/2; xx=linspace(a,b,500);
N=100; h=(b-a)/N; % N is the number of the iterations in the IVP
D=1; % the initial derivative to set properly
B=2*exp(1/2); % the value at the right hand endpoint

x(1)=a;

f=@(x,y)[y(2);y(1)+4*exp(x)];

y=[A;D];    

for i=1:N
    x(i+1)=x(i)+h;
    y(:,i+1)=y(:,i)+h*(1/2*f(x(i),y(:,i))+1/2*f(x(i)+h,y(:,i)+h*f(x(i),y(:,i))));
end

plot(x,y(1,:),'bo-',xx,yexact(xx),'r-')

%% Problem 76

yexact=@(x)(2*x+1).*exp(x);
a=0; A=1; b=1/2; xx=linspace(a,b,500);
N=3; h=(b-a)/(N+1); % N is the number of the inner points
x=[a+h:h:b-h]';

B=2*exp(1/2); % the value at the right hand endpoint

y=(gallery('tridiag',N,1,-2,1)/h^2-speye(N))\(4*exp(x)-[A/h^2,zeros(1,N-2),B/h^2]');
plot([a,x',b],[A,y',B],'bo-',xx,yexact(xx),'r-')

error=(norm(yexact(x)-y,inf))