%% Weeks 12-14
% Matlab codes for weeks 12-14

%% Problem 1

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 2

format long
n=4;
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 3

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 4

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=2 % 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 5

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

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 6

yexact=@(x)4*exp(x.^3/3)-1;
f=@(x,y)x.^2.*(1+y);
x0=0, y0=3,
xmax=1, N=8, 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 7

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 8

% 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 9

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 10

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 11

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-')

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