NUMERICAL TECHNICS

differentials: euler vs heun vs rungekutta

clear all;clf;
% euler method: y(i+1) = y(i) + df(t(i),y(i)) * h
tSample = 0.1;
tEnd = 10;

y = 2;
for t = 0:tSample:tEnd
    plot(t, y,'*r');
    hold on;
    y = y + (-2*y+4*exp(-t)) * tSample; 
end

% heun method: y0(i+1) = y(i) + df(t(i),y(i))*h and y(i+1) = y(i) + df(t(i),y(i)) + df(t(i+1),y0(i+1)) * h 
yh = 2;
for t = 0:tSample:tEnd
   	plot(t, yh,'*b');
  	hold on;
    
    y = yh;
    y0 = y + (-2*y+4*exp(-t)) * tSample; 
    yh = y + 0.5 * (((-2 * y + 4 * exp(-t)) + (-2 * y0 + 4 * exp(-t + tSample)))) * tSample; 
end

% 4th order RungeKutta:
x = [0:0.25:1];
f = @(xt,y) (1+2*xt)*sqrt(y);
[xt, y] = ode45(f, x, 1);                    % ode45 = 4de order RungeKutta
plot(xt, y)


% Example RungeKutta to how to use it to find max Y value:
clear all;clf;

% Exercise on differential with kutta runga:
g = 9.81;
R = 6370000;
t = [0:1:300];
v0 = 1400;

df = @(t,x) ((-g*R^2) / ((R+x)^2)) * t + v0;
[t,x] = ode45(df, t, v0);

yMax = 0;                                           % to find the top value
for i=1:length([t,x])
    if(x(i) > yMax)
        xPos = i;
        yMax = x(i);
    end
end
plot(t,x, xPos, yMax, '*r');

integrals: trapeziumrule vs simpsonrule

clear all;clf;

% Trapazium rule
x = [1 2 3.25 4.5 6 7 8 8.5 9.3 10];
y = [5 6 5.5 7 8.5 8 6 7 7 5];

t = linspace(0,10,1000);

p = polyfit(x,y,3);
b = polyval(p,t);

plot(x,y,'*r',t,b);

resultDist = polyval(polyint(p),x(length(x))) - polyval(polyint(p),x(1));

% simpsons rule

v =@(t) 1800*log(160000./(160000-2500*t)) - 9.81 * t;
a = 0;
b = 30;
result = ((b - a)/6) * (v(a) + 4*v(a+b/2) + v(b));

interpolate with polynomes vs matlab function

clear all;clf;
% interpolating with polyfit (5order polynome) and polyval:
x=[200,250,300,350,400,450];                % x is here the temp in Kelvin
y=[1.708,1.367,1.139,0.967,0.854,0.759];    % y is here de kg/m^3

xq = 330;                                   % we want to know value at 330° K

p = polyfit(x, y, 5);                       % interpolate with 5th order
b = polyval(p, xq);                         % polynome for 330° K
plot(x, y, xq, b, '*r');                            % plot it at x = 330°K

clear all;clf;
% matlab interpolate function:
x = [0,8,16,24,32,40];                          % x are temps in °C
y = [14.621,11.843,9.870,8.418,7.305,6.413];    % y are concentrations in mg/L

value = 27;                                    	% we want to estimate the concentration at 27°C

interLinear = interp1(x,y,value,'linear');
interSpline = interp1(x,y,value,'spline');

plot(x,y, value,interLinear,'*g', value, interSpline,'*r');

least squares method to plot data.

clear all;
% least squares method:
% input input x and output y.
x = [0 2 4 6 9 11 14 17 20];
y = [4 5 6 5 8 7 6 9 12];

a1 = (length(x) * sum(x.*y) - (sum(x)*sum(y))) / ((length(x) * sum(x.^2) - (sum(x)^2)));
a0 = mean(y) - a1 * mean(x);
    
yr = a0 + a1 * x;
plot (x, yr, x, y, '*r');

r = (length(x) * sum(x.*y) - sum(x) * sum(y)) / (sqrt(length(x) * sum(x.^2) - sum(x)^2) * sqrt(length(x) * sum(y.^2) - sum(y).^2)); 
r2 = r^2;                                   % check the correlation value of x and y.

Matrices: LU decomposition vs Gauss Seidel

clear all;
% LU Decomposition method:
A = [1 1 6; 1 5 -1; 4 2 -2];
b = [8; 5; 4];

[L,U] = lu(A);                      % matlab function is: x = A\b
x = U\(L\b);

% Gaussseidel method:
% 9x + 3y + z = 13; 6x + 8z = 2; 2x + 5y - z = 6
A = [9 3 1; 2 5 -1; 6 0 8];
B = [13; 6; 2];

iterMax = 1000;
for i=1:length(A)
    xold(i) = B(i)/A(1,i);
end
for iter=0:iterMax
    x(1) = (B(1) - A(1,2) * xold(2) - A(1,3) * xold(3)) / A(1,1);
    x(2) = (B(2) - A(2,1) * xold(1) - A(2,3) * xold(3)) / A(2,2);
    x(3) = (B(3) - A(3,1) * xold(1) - A(3,2) * xold(2)) / A(3,3);
    for i=1:length(A)
        xold(i) = x(i);
    end
end
answerGaussSeidel = x;

roots: fzero vs bisection vs RegulaFalsi vs NewtonRaphson vs Secant

clear all;

% Input formula:
f = @(x) 7*sin(x)*exp(-x)-1 
df = @(x) 7*exp(x)*(cos(x)-sin(x));

% matlabfunction:
matlabFunction = fzero(f,2);                    % random value to choose

% bisectie method:
iterMax = 1000;
xu = 2;
xl = 1;

if (f(xl) * f(xu) < 0)
    for iter=0:iterMax 
        xr = 0.5*(xl+xu);
        if (f(xr)*f(xl) < 0)
            xu = xr;
        else 
          	xl = xr;
        end  
        bisection = xr;
   end
end


% regula falsi method:
iterMax = 1000;
xu = 2;
xl = 1;

if ( f(xl)*f(xu) < 0 )
    for iter=0:iterMax
        xr = xu - ((f(xu)*(xl-xu))/(f(xl)-f(xu)));
        if ( f(xr)*f(xl) < 0 )
          	xu = xr;
        else
        	xl = xr;
        end 
        regulaFalsi = xr;
   end
end  


% newton raphson method:
iterMax = 1000;
x = 0;                                         	% start with x = 0;

x0 = 1                                        	% random value to choose
for iter=0:iterMax
    x0 = -f(x0)/df(x0) + x0;
    if ((x0-x)/x0 >= 0)
        x = x0;
    end
end
newtonRaphson = x0;

% secant method:
iterMax = 1000;
x = 1                                       	% start with x = 0;

x0 = 0.5;                                       % select fake random start value
x1 = 0.4;                                      	% select fake random start value
for iter=0:iterMax
    xk = x1 - f(x1) * ((x1-x0)/f(x1)-f(x0))
 	if (f(xk) ~= 0)                         
       	x0 = x1;
        x1 = xk;
    end
end  
secant = x1;