real

GaussPivot.m

@@ -0,0 +1,22 @@
+function x = GaussPivot(A,b)
+[m,n]=size(A);
+if m~=n, error('Matrix A must be square'); end
+nb=n+1;
+Aug=[A b];
+for k = 1:n-1
+ [big,i]=max(abs(Aug(k:n,k)));
+ ipr=i+k-1;
+ if ipr~=k
+ Aug([k,ipr],:)=Aug([ipr,k],:);
+ end
+ for i = k+1:n
+  factor=Aug(i,k)/Aug(k,k);
+  Aug(i,k:nb)=Aug(i,k:nb)-factor*Aug(k,k:nb);
+ end
+end
+x=zeros(n,1);
+x(n)=Aug(n,nb)/Aug(n,n);
+for i = n-1:-1:1
+ x(i)=(Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);
+end
+end

LU_pivot.m

@@ -0,0 +1,18 @@
+function [L, U, P] = LU_pivot(A)
+[n, ~] = size(A); L = eye(n); P = eye(n); U = A;
+for j = 1:n
+ [~, m] = max(abs(U(j:n, j)));
+ m = m + j - 1;
+ if m ~= j
+  U([m,j], :) = U([j,m], :);
+  P([m,j], :) = P([j,m], :);
+  if j >= 2
+   L([m,j], 1:j-1) = L([j,m], 1:j-1);
+  end
+ end
+ for i = j+1:n
+  L(i, j) = U(i, j) / U(j, j);
+  U(i, :) = U(i, :) - L(i, j)*U(j, :);
+ end
+end
+end

PVSSR.m

@@ -0,0 +1,37 @@
+function f = PVSSR(a, x, y)
+% PVSSR: Calculates the sum of squared residuals for the real gas equation.
+%   f = PVSSR(a, x, y)
+%   a: vector of parameters [A1, A2]
+%   x: independent variable (Volume P in this case)
+%   y: dependent variable (P in this case)
+%   f: Sum of Squared Residuals (SSR)
+
+% Constants from the problem
+R = 82.05; % mL.atm/gmol.K
+T = 303;   % K
+
+if length(x) ~= length(y)
+    fprintf('The two data must have the same length\n');
+    f = NaN; % Return NaN if dimensions don't match
+    return;
+end
+
+n = length(x); % The number of data
+
+yp_model = zeros(1, n);
+for j = 1:n
+    % The equation: P = RT/V * (1 + A1/V + A2/V^2)
+    % Here, x(j) is Volume (V)
+    % a(1) is A1
+    % a(2) is A2
+    yp_model(j) = R * T * (1 + a(1)/x(j) + a(2)/x(j)^2) / x(j);
+end
+
+% Calculate the sum of squared residuals
+% f = (y - yp_model) * (y - yp_model)';
+% This is the sum of squares. Ensure y and yp_model are row vectors.
+y = reshape(y, 1, []); % Ensure y is a row vector
+yp_model = reshape(yp_model, 1, []); % Ensure yp_model is a row vector
+
+f = sum((y - yp_model).^2); % Sum of squared residuals
+end

bisection.m

@@ -0,0 +1,38 @@
+function [root,fx,ea,iter]=bisection(func,xl,xu,es,maxit)
+% bisect: root location zeroes
+% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
+% uses bisection method to find the root of func
+% input:
+% func = name of function
+% xl, xu = lower and upper guesses
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by func
+% output:
+% root = real root
+% fx = function value at root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+test = func(xl)*func(xu);
+if test>0,error('no sign change'),end
+if nargin<4||isempty(es), es=0.0001;end
+if nargin<5||isempty(maxit), maxit=50;end
+iter = 0; xr = xl; ea = 100;
+while (1)
+xrold = xr;
+xr = (xl + xu)/2;
+iter = iter + 1;
+if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
+test = func(xl)*func(xr);
+if test < 0
+xu = xr;
+elseif test > 0
+xl = xr;
+else
+ea = 0;
+end
+if ea <= es || iter >= maxit,break,end
+end
+root = xr;
+fx = func(xr);

lab2/lab2.m

@@ -0,0 +1,206 @@
+%% Labwork 2 Solutions
+
+%%% III. Polynomials
+
+%% Exercise 1: Arithmetic Operations
+
+% a) 12^2 - 13^2
+a = 12^2 - 13^2
+
+% b) 24 * 5 / 1.5
+b = 24 * 5 / 1.5
+
+% c) 7! (7 factorial)
+c = factorial(7)
+
+% d) 4^2 + 3^3 - 124
+d = 4^2 + 3^3 - 124
+
+% e) 6/5 * (1/2 * 1/4 - 1)
+e = (6/5) * ((1/2) * (1/4) - 1)
+
+%% Exercise 2: fzero
+
+% Given f(x) = x^2 - 9
+
+% a) Find the positive root using fzero
+f = @(x) x^2 - 9;  % Define the function
+positive_root = fzero(f, 2)  % Initial guess of 2 (can be any positive number)
+
+% b) Initial guess of zero
+root_zero_guess = fzero(f, 0) % It might give a warning, and the result can vary slightly due to how fzero works
+
+% c) Display all iterations
+options = optimset('Display', 'iter');
+[root_iterations, fval_iterations, exitflag_iterations, output_iterations] = fzero(f, 2, options)
+
+% d) Set tolerance: tolx = 10^-3
+options_tol = optimset('TolX', 10^-3);
+root_tolerance = fzero(f, 2, options_tol)
+
+% e) Verify with roots
+coefficients = [1 0 -9]; % Coefficients of the quadratic equation
+all_roots = roots(coefficients)
+
+
+%% Exercise 3: plot
+
+% Given v(t) = sqrt(gm/cd) * tanh(sqrt(g*cd/m) * t)
+% g = 9.81, m = 68.1, cd = 0.25
+
+g = 9.81;
+m = 68.1;
+cd = 0.25;
+
+t = 0:2:20; % Time vector from 0 to 20 seconds in steps of 2
+v = sqrt((g*m)/cd) .* tanh(sqrt((g*cd)/m) .* t);  % Calculate velocity
+
+% a) Plot (t, v)
+figure; % Create a new figure
+plot(t, v, 'b-o'); % Plot v against t ('b-o' for blue line with circles)
+
+% b) Customize the graph
+title('Plot of v versus t');
+xlabel('Values of t');
+ylabel('Values of v');
+grid on;
+
+
+%%% IV. System of Non-linear Equations
+
+%% Exercise 1: Newton-Raphson
+
+% a) Use ezplot to derive good guesses
+%   (This part is visual, so I'll show the code to generate the plot for you to estimate)
+figure;
+ezplot(@(x1, x2) x1^2 + x1*x2 - 10, [-6 6]); % Plot the first equation
+hold on;
+ezplot(@(x1, x2) x2 + 3*x1*x2^2 - 57, [-6 6]); % Plot the second equation
+hold off;
+title('Plot of equations for initial guesses');
+xlabel('x1');
+ylabel('x2');
+grid on;
+%   Inspect the plot to get approximate intersection points (guesses).
+
+% b) Explain the given code
+%   (Explanation embedded as comments in the code below)
+
+clc  % Clear command window
+clear all  % Clear variables
+close all  % Close all figures
+
+% Define the initial guess as a column vector
+x = [1; 1];  
+
+N = 50;  % Maximum number of iterations
+
+for k = 1:N  % Loop for iterations
+    % Define the system of equations
+    F = [x(1)^2 + x(1)*x(2) - 10; 
+         x(2) + 3*x(1)*x(2)^2 - 57];  
+    
+    % Calculate the Jacobian matrix
+    A = [2*x(1) + x(2), x(1); 
+         3*x(2)^2, 1 + 6*x(1)*x(2)];  
+    
+    % Solve for the change in x (dx)
+    dx = A \ F;  
+    
+    % Update x
+    x = x - dx;  
+    
+    % Increment the iteration counter (though it's redundant here, as the loop does it)
+    %k = k + 1; 
+end
+
+k  % Display the final iteration number
+x'  % Display the solution as a row vector
+F  % Display the function values at the solution (should be close to zero)
+
+
+% c) Create root2d.m and use fsolve
+%
+%   Create a new MATLAB function file named root2d.m with the following content:
+%   ------------------------------------------------------------------
+%   function F = root2d(x)
+%       F(1) = x(1)^2 + x(1)*x(2) - 10;
+%       F(2) = x(2) + 3*x(1)*x(2)^2 - 57;
+%   end
+%   ------------------------------------------------------------------
+%
+%   Then, in a separate script:
+
+% Call fsolve with the function and an initial guess
+f_system = @root2d;
+x0 = [0; 0];  % Initial guess
+[solution_fsolve, fval_fsolve] = fsolve(f_system, x0)
+
+%   Explanation:
+%   - root2d.m defines the system of equations as a function that takes a vector x 
+%     and returns a vector F representing the equations.
+%   - fsolve finds the solution to the system of nonlinear equations.  It takes the
+%     function handle (@root2d) and an initial guess (x0) as input.
+
+% d) Use newtmult.m
+%   (Assuming newtmult.m is provided as in the lab, place it in the same directory)
+
+% Define the function (same as in root2d.m)
+func = @(x) [x(1)^2 + x(1)*x(2) - 10; x(2) + 3*x(1)*x(2)^2 - 57];
+x0 = [0; 0];  % Initial guess
+es = 0.0001; % Desired relative error
+maxit = 50;  % Maximum iterations
+
+[x_newtmult, f_newtmult, ea_newtmult, iter_newtmult] = newtmult(func, x0, es, maxit)
+
+%   Explanation of newtmult.m (from the lab):
+%   - It implements the Newton-Raphson method for systems of nonlinear equations.
+%   - Inputs: function handle (func), initial guess (x0), desired relative error (es), 
+%     maximum iterations (maxit), and any additional parameters for the function.
+%   - Outputs: solution (x), function values at the solution (f), approximate relative error (ea), 
+%     and number of iterations (iter).
+%   - The code iteratively calculates the Jacobian, solves a linear system, and updates the 
+%     solution until convergence or the maximum number of iterations is reached.
+
+
+%%% V. Common
+
+%% Exercise 1: Calculate the value of x
+
+% a) log(x) = 770  (Assuming base 10 logarithm)
+x_a = 10^770
+
+% b) e^x = 15
+x_b = log(15)
+
+% c) sin(x) = 0.12 (in radians)
+x_c = asin(0.12)
+
+% d) sin(x + pi/5) = 0.25 (in radians)
+x_d = asin(0.25) - pi/5
+
+% e) tan(x + 2*pi/3) = 3 (in radians)
+x_e = atan(3) - 2*pi/3
+
+%% Exercise 2: Create a vector x
+
+% Create the row vector
+x_vector = [2 3 1 5 11]
+
+% a) Find length and size
+length_x = length(x_vector)
+size_x = size(x_vector)
+%   - length() returns the number of elements in the longest dimension.
+%   - size() returns a row vector with the number of rows and columns.
+
+% b) Find the transpose
+x_transpose = x_vector'
+
+% c) Calculate the product of x and its transpose
+product_x_xt = x_vector * x_transpose
+
+
+%% Exercise 3: Build the matrix B
+
+% Build the matrix B
+B = [1:7; 9:-2:-3; 4.^(1:7)]

lab3_old/GaussNoPivot.m

@@ -0,0 +1,18 @@
+% Modified Gauss Elimination without pivoting
+function x = GaussNoPivot(A,b)
+[m,n]=size(A);
+if m~=n, error('Matrix A must be square'); end
+nb=n+1;
+Aug=[A b];
+for k = 1:n-1
+ for i = k+1:n
+  factor=Aug(i,k)/Aug(k,k);
+  Aug(i,k:nb)=Aug(i,k:nb)-factor*Aug(k,k:nb);
+ end
+end
+x=zeros(n,1);
+x(n)=Aug(n,nb)/Aug(n,n);
+for i = n-1:-1:1
+ x(i)=(Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);
+end
+end

lab3_old/mylup.m

@@ -0,0 +1,18 @@
+function [L, U, P] = mylup(A)
+[n,~] = size(A);
+C = [A (1:n)'];
+for j = 1:n-1
+ [~, i] = max(abs(C(j:n, j)));
+ i = i + j - 1;
+ if i ~= j
+  C([i j], :) = C([j i], :);
+ end
+ for k = j+1:n
+  C(k, j) = C(k, j)/C(j, j);
+  C(k, j+1:n) = C(k, j+1:n) - C(k, j)*C(j, j+1:n);
+ end
+end
+PP = C(:, end); C = C(:, 1:n);
+L = tril(C, -1) + eye(n); U = triu(C);
+P = eye(n); P = P(PP, :);
+end

lab4/createFit.m

@@ -0,0 +1,35 @@
+function [fitresult, gof] = createFit(mass_kg, metabolism_watts)
+%CREATEFIT(MASS_KG,METABOLISM_WATTS)
+%  Create a fit.
+%
+%  Data for 'untitled fit 1' fit:
+%      X Input: mass_kg
+%      Y Output: metabolism_watts
+%  Output:
+%      fitresult : a fit object representing the fit.
+%      gof : structure with goodness-of fit info.
+%
+%  See also FIT, CFIT, SFIT.
+
+%  Auto-generated by MATLAB on 29-May-2025 01:56:21
+
+
+%% Fit: 'untitled fit 1'.
+[xData, yData] = prepareCurveData( mass_kg, metabolism_watts );
+
+% Set up fittype and options.
+ft = fittype( 'poly2' );
+
+% Fit model to data.
+[fitresult, gof] = fit( xData, yData, ft, 'Normalize', 'on' );
+
+% Plot fit with data.
+figure( 'Name', 'untitled fit 1' );
+h = plot( fitresult, xData, yData );
+legend( h, 'metabolism_watts vs. mass_kg', 'untitled fit 1', 'Location', 'NorthEast', 'Interpreter', 'none' );
+% Label axes
+xlabel( 'mass_kg', 'Interpreter', 'none' );
+ylabel( 'metabolism_watts', 'Interpreter', 'none' );
+grid on
+
+

lab4/l4ame.m

@@ -0,0 +1,13 @@
+real = [40 1450
+63 3817
+62 3708
+68 4300
+64 3950
+45 1896
+50 2339
+66 4200
+67 4312]
+
+f = polyfit(real(:, 1), real(:, 2), 1)
+
+disp(polyval(f, 72))

lab4/lab4a.m

@@ -0,0 +1,386 @@
+%% Labwork 4: Curve Fitting Problems
+
+%%% I. Curve Fitting
+
+%% Exercise 1: Data Fitting with Various Functions
+
+% Data from the table
+weight = [40; 63; 62; 68; 64; 45; 50; 66; 67]; % kg
+jhi = [1450; 3817; 3708; 4300; 3950; 1896; 2339; 4200; 4312]; % JHI number
+
+% Weight of Van Lam to predict
+van_lam_weight = 72;
+
+% 1. Use function "fit()" in Matlab. Plot results
+
+% Define the type of fit. A linear fit is often a good starting point for such data.
+ft = fittype('poly1'); % Linear polynomial fit (y = ax + b)
+
+% Fit the data
+[fit_obj_fit, gof_fit] = fit(weight, jhi, ft);
+
+% Predict JHI for Van Lam
+jhi_van_lam_fit = fit_obj_fit(van_lam_weight);
+
+% Plot results
+figure;
+plot(fit_obj_fit, weight, jhi); % Plot the fit and the data points
+hold on;
+plot(van_lam_weight, jhi_van_lam_fit, 'rx', 'MarkerSize', 10, 'LineWidth', 2); % Mark Van Lam's prediction
+xlabel('Weight (kg)');
+ylabel('JHI Number');
+title('JHI vs. Weight using fit()');
+legend('Data', 'Linear Fit', 'Van Lam''s Prediction');
+grid on;
+hold off;
+
+fprintf('1. Predicted JHI for Van Lam using fit(): %.2f\n', jhi_van_lam_fit);
+
+% 2. Use function "polyfit()" in Matlab. Plot results
+
+% Fit a first-degree polynomial (linear fit)
+degree = 1;
+coefficients_polyfit = polyfit(weight, jhi, degree);
+
+% Create a polynomial function from coefficients
+poly_fit_func = @(x) coefficients_polyfit(1) * x + coefficients_polyfit(2);
+
+% Predict JHI for Van Lam
+jhi_van_lam_polyfit = poly_fit_func(van_lam_weight);
+
+% Plot results
+figure;
+plot(weight, jhi, 'o'); % Plot data points
+hold on;
+x_plot = linspace(min(weight), max(weight) + 10, 100); % Extend plot range slightly
+plot(x_plot, poly_fit_func(x_plot), 'b-'); % Plot the fitted line
+plot(van_lam_weight, jhi_van_lam_polyfit, 'rx', 'MarkerSize', 10, 'LineWidth', 2); % Mark Van Lam's prediction
+xlabel('Weight (kg)');
+ylabel('JHI Number');
+title('JHI vs. Weight using polyfit()');
+legend('Data', 'Linear Fit', 'Van Lam''s Prediction');
+grid on;
+hold off;
+
+fprintf('2. Predicted JHI for Van Lam using polyfit(): %.2f\n', jhi_van_lam_polyfit);
+
+
+% 3. Use function "lsqcurvefit()" in Matlab. Plot results
+
+% Define the model function for lsqcurvefit (linear model: y = p(1)*x + p(2))
+model_func = @(p, x) p(1) * x + p(2);
+
+% Initial guess for parameters [slope, intercept]
+p0 = [10; 0]; % Arbitrary initial guess
+
+% Perform the least-squares curve fit
+[parameters_lsq, resnorm_lsq] = lsqcurvefit(model_func, p0, weight, jhi);
+
+% Predict JHI for Van Lam
+jhi_van_lam_lsqcurvefit = model_func(parameters_lsq, van_lam_weight);
+
+% Plot results
+figure;
+plot(weight, jhi, 'o'); % Plot data points
+hold on;
+x_plot = linspace(min(weight), max(weight) + 10, 100); % Extend plot range
+plot(x_plot, model_func(parameters_lsq, x_plot), 'g-'); % Plot the fitted line
+plot(van_lam_weight, jhi_van_lam_lsqcurvefit, 'rx', 'MarkerSize', 10, 'LineWidth', 2); % Mark Van Lam's prediction
+xlabel('Weight (kg)');
+ylabel('JHI Number');
+title('JHI vs. Weight using lsqcurvefit()');
+legend('Data', 'Linear Fit', 'Van Lam''s Prediction');
+grid on;
+hold off;
+
+fprintf('3. Predicted JHI for Van Lam using lsqcurvefit(): %.2f\n', jhi_van_lam_lsqcurvefit);
+
+
+% 4. Give some comments for the differences among these 3 methods.
+fprintf('\n4. Comments on the differences among fit(), polyfit(), and lsqcurvefit():\n');
+fprintf('   - fit(): A high-level function from the Curve Fitting Toolbox. It provides a more generalized framework\n');
+fprintf('            for curve fitting, allowing for various predefined models (polynomials, exponentials, etc.)\n');
+fprintf('            and custom equations. It returns a "fittype" object, which is convenient for plotting and prediction.\n');
+fprintf('            It also provides goodness-of-fit statistics (gof).\n');
+fprintf('   - polyfit(): A basic MATLAB function specifically designed for polynomial regression. It returns the\n');
+fprintf('                polynomial coefficients. It is straightforward for polynomial fits but less flexible for\n');
+fprintf('                non-polynomial or custom models compared to `fit()` or `lsqcurvefit()`.\n');
+fprintf('   - lsqcurvefit(): A function from the Optimization Toolbox designed for solving non-linear least-squares\n');
+fprintf('                    problems. It requires you to define the model function explicitly and provide an initial guess\n');
+fprintf('                    for the parameters. It is highly flexible for custom non-linear models and allows for bounds\n');
+fprintf('                    and options for the optimization process. It is more fundamental for custom curve fitting.\n');
+fprintf('   - For this linear problem, all three methods yield very similar results as they are all performing\n');
+fprintf('     least-squares linear regression under the hood. The choice depends on the complexity of the model\n');
+fprintf('     and desired level of control/toolboxes available.\n');
+
+
+%% Exercise 2: Non-linear Curve Fitting (lsqnonlin)
+
+% Given decaying equation: y = exp(-t.x) + 0.8
+% Note: 't' here is the parameter to find, 'x' is the independent variable.
+% The problem uses 't' as a parameter and 'x' as independent variable, which can be confusing.
+% Let's use 'p' for the parameter (decaying rate) to avoid confusion with time 't'.
+
+% Recorded observations
+x_obs = [0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1.0];
+y_obs = [1.6781; 1.5711; 1.4771; 1.3945; 1.3220; 1.2584; 1.2025; 1.1535; 1.1104; 1.0725];
+
+% Define the objective function for lsqnonlin: residual = F(x, p) - y_observed
+% fun = @(p) (exp(-p(1) * x_obs) + 0.8) - y_obs;
+% The above defines the residual. lsqnonlin minimizes the sum of squares of these residuals.
+% Let's redefine 'fun' to represent the model and use lsqnonlin's direct form
+% to pass the data.
+
+% The problem statement is slightly ambiguous with 't' and 'x'.
+% Assuming 't' in the equation y = exp(-t.x)+0.8 is the *parameter* to be found.
+% Let's rename the parameter to `decay_rate` for clarity.
+
+% Objective function for lsqnonlin (returns the residuals)
+% The input `p` is the parameter (decaying rate)
+fun = @(p) (exp(-p(1) * x_obs) + 0.8) - y_obs;
+
+% a. Find the decaying rate t by function "lsqnonlin()" in Matlab. Plot results.
+
+% Initial guess for the decaying rate (t in the problem, p(1) here)
+x0_lsqnonlin = 1; % Initial guess for 't' (decaying rate)
+
+% Perform the non-linear least squares optimization
+[decay_rate_found, resnorm_lsqnonlin, residual_lsqnonlin, exitflag_lsqnonlin, output_lsqnonlin] = lsqnonlin(fun, x0_lsqnonlin);
+
+% Reconstruct the fitted curve using the found decay rate
+y_fitted = exp(-decay_rate_found * x_obs) + 0.8;
+
+% Plot results
+figure;
+plot(x_obs, y_obs, 'o', 'DisplayName', 'Observed Data'); % Plot observed data
+hold on;
+plot(x_obs, y_fitted, 'r-', 'DisplayName', ['Fitted Curve (t = ' num2str(decay_rate_found) ')']); % Plot fitted curve
+xlabel('x');
+ylabel('y');
+title('Kagawa''s Particle Decaying Equation');
+legend('show');
+grid on;
+hold off;
+
+fprintf('\nExercise 2a: Found decaying rate (t) using lsqnonlin(): %.4f\n', decay_rate_found);
+
+
+% b. Try to use lb, ub, options in: x = lsqnonlin(fun,x0,lb,ub,options)
+
+% Define lower bound (lb) and upper bound (ub) for the parameter
+% Assuming the decaying rate 't' should be positive.
+lb = 0;   % Lower bound for 't'
+ub = 5; % Upper bound for 't' (arbitrarily chosen, based on expected range)
+
+% Define optimization options
+options_lsqnonlin = optimoptions('lsqnonlin', ...
+                                 'Display', 'iter', ... % Display iteration progress
+                                 'Algorithm', 'levenberg-marquardt', ... % Specify algorithm
+                                 'TolX', 1e-6); % Tolerance on x
+
+[decay_rate_bounded, resnorm_bounded] = lsqnonlin(fun, x0_lsqnonlin, lb, ub, options_lsqnonlin);
+
+fprintf('\nExercise 2b: Found decaying rate (t) with bounds and options: %.4f\n', decay_rate_bounded);
+
+
+%% Exercise 3: Interpolation (interp1)
+
+% Given table of observation
+x_obs_interp = [0.00000; 0.78540; 1.57080; 2.35619; 3.14159; 3.92699; 4.71239; 5.49779; 6.28319];
+f_obs_interp = [0.00000; 0.70711; 1.00000; 0.70711; 0.00000; -0.70711; -1.00000; -0.70711; 0.00000];
+
+% Points to interpolate
+xq1 = 3.00000;
+xq2 = 4.50000;
+
+% a. Calculate f(3.00000) and f(4.50000) by function "interp1()" in Matlab. Plot results.
+
+% Default method for interp1 is 'linear'
+f_xq1_linear = interp1(x_obs_interp, f_obs_interp, xq1);
+f_xq2_linear = interp1(x_obs_interp, f_obs_interp, xq2);
+
+% Plot results
+figure;
+plot(x_obs_interp, f_obs_interp, 'o', 'DisplayName', 'Observed Data'); % Original data
+hold on;
+x_interp_plot = linspace(min(x_obs_interp), max(x_obs_interp), 100);
+f_interp_linear = interp1(x_obs_interp, f_obs_interp, x_interp_plot, 'linear');
+plot(x_interp_plot, f_interp_linear, 'b-', 'DisplayName', 'Linear Interpolation'); % Linear interpolation curve
+plot(xq1, f_xq1_linear, 'gx', 'MarkerSize', 10, 'LineWidth', 2, 'DisplayName', 'f(3.0) (Linear)'); % Interpolated points
+plot(xq2, f_xq2_linear, 'mx', 'MarkerSize', 10, 'LineWidth', 2, 'DisplayName', 'f(4.5) (Linear)');
+xlabel('x');
+ylabel('f(x)');
+title('Interpolation using interp1() (Linear)');
+legend('show');
+grid on;
+hold off;
+
+fprintf('\nExercise 3a: f(3.00000) using linear interpolation: %.5f\n', f_xq1_linear);
+fprintf('Exercise 3a: f(4.50000) using linear interpolation: %.5f\n', f_xq2_linear);
+
+
+% b. Try to use another method ('nearest', 'cubic', or 'spline'.) with vq = interp1(x,v,xq,method)
+
+% Interpolate using 'nearest'
+f_xq1_nearest = interp1(x_obs_interp, f_obs_interp, xq1, 'nearest');
+f_xq2_nearest = interp1(x_obs_interp, f_obs_interp, xq2, 'nearest');
+fprintf('\nExercise 3b: f(3.00000) using nearest interpolation: %.5f\n', f_xq1_nearest);
+fprintf('Exercise 3b: f(4.50000) using nearest interpolation: %.5f\n', f_xq2_nearest);
+
+% Interpolate using 'cubic'
+f_xq1_cubic = interp1(x_obs_interp, f_obs_interp, xq1, 'cubic');
+f_xq2_cubic = interp1(x_obs_interp, f_obs_interp, xq2, 'cubic');
+fprintf('Exercise 3b: f(3.00000) using cubic interpolation: %.5f\n', f_xq1_cubic);
+fprintf('Exercise 3b: f(4.50000) using cubic interpolation: %.5f\n', f_xq2_cubic);
+
+% Interpolate using 'spline'
+f_xq1_spline = interp1(x_obs_interp, f_obs_interp, xq1, 'spline');
+f_xq2_spline = interp1(x_obs_interp, f_obs_interp, xq2, 'spline');
+fprintf('Exercise 3b: f(3.00000) using spline interpolation: %.5f\n', f_xq1_spline);
+fprintf('Exercise 3b: f(4.50000) using spline interpolation: %.5f\n', f_xq2_spline);
+
+% Plot all interpolation methods for comparison
+figure;
+plot(x_obs_interp, f_obs_interp, 'o', 'DisplayName', 'Observed Data');
+hold on;
+plot(x_interp_plot, interp1(x_obs_interp, f_obs_interp, x_interp_plot, 'linear'), 'b-', 'DisplayName', 'Linear');
+plot(x_interp_plot, interp1(x_obs_interp, f_obs_interp, x_interp_plot, 'nearest'), 'g--', 'DisplayName', 'Nearest');
+plot(x_interp_plot, interp1(x_obs_interp, f_obs_interp, x_interp_plot, 'cubic'), 'r:', 'DisplayName', 'Cubic');
+plot(x_interp_plot, interp1(x_obs_interp, f_obs_interp, x_interp_plot, 'spline'), 'm-.', 'DisplayName', 'Spline');
+
+plot(xq1, f_xq1_linear, 'b^', 'MarkerSize', 8, 'DisplayName', 'f(3.0) Linear');
+plot(xq1, f_xq1_nearest, 'gv', 'MarkerSize', 8, 'DisplayName', 'f(3.0) Nearest');
+plot(xq1, f_xq1_cubic, 'r>', 'MarkerSize', 8, 'DisplayName', 'f(3.0) Cubic');
+plot(xq1, f_xq1_spline, 'm<', 'MarkerSize', 8, 'DisplayName', 'f(3.0) Spline');
+
+plot(xq2, f_xq2_linear, 'b*', 'MarkerSize', 8, 'DisplayName', 'f(4.5) Linear');
+plot(xq2, f_xq2_nearest, 'gs', 'MarkerSize', 8, 'DisplayName', 'f(4.5) Nearest');
+plot(xq2, f_xq2_cubic, 'rd', 'MarkerSize', 8, 'DisplayName', 'f(4.5) Cubic');
+plot(xq2, f_xq2_spline, 'mh', 'MarkerSize', 8, 'DisplayName', 'f(4.5) Spline');
+
+xlabel('x');
+ylabel('f(x)');
+title('Comparison of Interpolation Methods');
+legend('show');
+grid on;
+hold off;
+
+
+%% Exercise 4: 2D Interpolation (interp2) and Surface Plotting
+
+% Given x runs from -3 to 3 with step of 1 and so does y.
+x_grid_base = -3:1:3;
+y_grid_base = -3:1:3;
+
+% 1. Draw a grid of coordinates of (x,y) (hint: meshgrid)
+[X, Y] = meshgrid(x_grid_base, y_grid_base);
+
+figure;
+plot(X, Y, 'b.'); % Plot grid points as blue dots
+hold on;
+plot(X', Y', 'b.'); % Also plot transposed to get vertical lines of dots
+xlabel('x');
+ylabel('y');
+title('Grid of (x,y) Coordinates');
+grid on;
+axis equal; % Equal scaling for x and y axes
+hold off;
+
+
+% 2. Build a random function z having many local extremes that takes only 2 arguments x and y. (hint: peaks)
+Z = peaks(X, Y);
+
+figure;
+mesh(X, Y, Z); % Plot as a mesh surface
+xlabel('x');
+ylabel('y');
+zlabel('z');
+title('Function z = peaks(X, Y)');
+colorbar; % Add a color bar
+grid on;
+
+
+% 3. Compute z(1.2, 2.8) and z(2.5, 2.5) by function "interp2()"
+
+% Define query points
+xq = [1.2, 2.5];
+yq = [2.8, 2.5];
+
+% Perform 2D interpolation (default method is 'linear')
+zq = interp2(X, Y, Z, xq, yq);
+
+fprintf('\nExercise 4: Interpolated Z values:\n');
+fprintf('z(1.2, 2.8) = %.4f\n', zq(1));
+fprintf('z(2.5, 2.5) = %.4f\n', zq(2));
+
+% Plot the interpolated points on the surface
+figure;
+mesh(X, Y, Z); % Original surface
+hold on;
+plot3(xq(1), yq(1), zq(1), 'ro', 'MarkerSize', 10, 'LineWidth', 2, 'MarkerFaceColor', 'r'); % Plot first point
+plot3(xq(2), yq(2), zq(2), 'gx', 'MarkerSize', 10, 'LineWidth', 2, 'MarkerFaceColor', 'g'); % Plot second point
+xlabel('x');
+ylabel('y');
+zlabel('z');
+title('Function z = peaks(X, Y) with Interpolated Points');
+legend('Surface', 'z(1.2, 2.8)', 'z(2.5, 2.5)');
+colorbar;
+grid on;
+hold off;
+
+
+%%% II. Common
+
+%% Exercise 1. Find the eigenvalues and eigenvectors of the following matrix
+%  (2 1)
+%  (1 2)
+
+A_eig = [2 1; 1 2];
+
+% Compute eigenvalues (D) and eigenvectors (V)
+[V, D] = eig(A_eig);
+
+fprintf('\nExercise II.1: Eigenvalues and Eigenvectors of matrix A:\n');
+fprintf('Eigenvalues (diagonal of D):\n');
+disp(diag(D));
+fprintf('Eigenvectors (columns of V):\n');
+disp(V);
+
+
+%% Exercise 2. Plot data
+
+% Given data
+T_celsius = [5 20 30 50 55]; % Temperature (C)
+mu = [0.08 0.015 0.009 0.006 0.0055]; % Viscosity (mu)
+
+% Plot with '*' mode
+figure;
+plot(T_celsius, mu, '*');
+xlabel('Temperature (C)');
+ylabel('Viscosity (mu)');
+title('Viscosity vs. Temperature (Marker: *)');
+grid on;
+
+% Plot with 'o' mode
+figure;
+plot(T_celsius, mu, 'o');
+xlabel('Temperature (C)');
+ylabel('Viscosity (mu)');
+title('Viscosity vs. Temperature (Marker: o)');
+grid on;
+
+% Plot with '.' mode
+figure;
+plot(T_celsius, mu, '.');
+xlabel('Temperature (C)');
+ylabel('Viscosity (mu)');
+title('Viscosity vs. Temperature (Marker: .)');
+grid on;
+
+% Plot with lines and markers for clarity (optional, but good practice)
+figure;
+plot(T_celsius, mu, 'b-o', 'LineWidth', 1.5, 'MarkerSize', 6, 'MarkerFaceColor', 'b');
+xlabel('Temperature (C)');
+ylabel('Viscosity (mu)');
+title('Viscosity vs. Temperature (Line and Markers)');
+grid on;

lab4/lab4b.m

@@ -0,0 +1,328 @@
+%% Labwork 4 (Part B): Regression and Statistics
+
+%%% 1. Linear Regression
+
+% Data
+mass_kg = [400; 70; 45; 2; 0.3; 0.16]; % Independent variable (x)
+metabolism_watts = [270; 82; 50; 4.8; 1.45; 0.97]; % Dependent variable (y)
+
+% Tiger's mass
+tiger_mass = 350; % kg
+
+%% (a) Write a script file which employs the following function file: linear_regression.m
+
+% Make sure linear_regression.m is in the current directory or on MATLAB's path.
+% Call the custom linear regression function
+[a0_custom, a1_custom, r2_custom] = linear_regression(mass_kg, metabolism_watts);
+
+% Predict metabolism rate for the 350-kg tiger
+predicted_metabolism_custom = a0_custom + a1_custom * tiger_mass;
+
+fprintf('--- 1. Linear Regression (Custom Function) ---\n');
+fprintf('  a0 (intercept) = %.4f\n', a0_custom);
+fprintf('  a1 (slope) = %.4f\n', a1_custom);
+fprintf('  r^2 (coefficient of determination) = %.4f\n', r2_custom);
+fprintf('  Predicted metabolism for 350kg tiger: %.4f watts\n', predicted_metabolism_custom);
+
+% The linear_regression function already plots, but we can enhance it if needed
+figure(1); % Bring the plot generated by linear_regression to front
+hold on;
+plot(tiger_mass, predicted_metabolism_custom, 'ro', 'MarkerSize', 8, 'LineWidth', 2, 'DisplayName', 'Tiger Prediction');
+legend('show');
+title('Linear Regression (Custom Function) with Tiger Prediction');
+hold off;
+
+
+%% (b) Use the Curve fitting Apps in MATLAB. Is there any difference in the two regression results?
+
+fprintf('\n--- 1. Linear Regression (Curve Fitting App) ---\n');
+fprintf('  To perform this:\n');
+fprintf('  1. Go to the "APPS" tab in the MATLAB toolstrip.\n');
+fprintf('  2. Find and click "Curve Fitting".\n');
+fprintf('  3. In the Curve Fitter app, go to "Data" -> "Select Data".\n');
+fprintf('  4. Select ''mass_kg'' for X data and ''metabolism_watts'' for Y data.\n');
+fprintf('  5. For "Fit type", select "Linear".\n');
+fprintf('  6. Observe the coefficients (p1 and p2, corresponding to a1 and a0) and R-squared.\n');
+fprintf('  7. You can then use the "Evaluate" panel to predict for 350kg.\n');
+fprintf('  Comparison: Results should be identical to the custom function if the custom function is correctly implemented\n');
+fprintf('  for standard linear regression, as both perform the same least-squares calculation.\n');
+
+
+%% (c) Use the functions polyfit in MATLAB to fit the points by a polynomial with degree n and polyval to estimate the metabolism rate of the tiger.
+% From what degree of the polynomial the estimated result becomes badly conditioned?
+% And from what degree of the polynomial the fitting curve is no longer unique?
+
+fprintf('\n--- 1. Linear Regression (polyfit/polyval) ---\n');
+
+% We'll try different degrees and observe behavior
+degrees_to_test = 1:min(length(mass_kg)-1, 5); % Test up to degree 5 or n-1
+
+for n_degree = degrees_to_test
+    fprintf('\n  --- Degree %d Polynomial Fit ---\n', n_degree);
+    
+    % Polyfit for coefficients
+    % polyfit can issue warnings for ill-conditioned polynomials at high degrees
+    warning('off', 'MATLAB:polyfit:RepeatedRoots'); % Suppress warnings about repeated roots (can happen for high degrees)
+    warning('off', 'MATLAB:polyfit:PolyNotUnique'); % Suppress warnings about non-unique poly
+    
+    try
+        coefficients_poly = polyfit(mass_kg, metabolism_watts, n_degree);
+        
+        % Predict metabolism for the tiger
+        predicted_metabolism_poly = polyval(coefficients_poly, tiger_mass);
+        
+        fprintf('  Coefficients for degree %d: [', n_degree);
+        fprintf(' %.4f', coefficients_poly);
+        fprintf(' ]\n');
+        fprintf('  Predicted metabolism for 350kg tiger (degree %d): %.4f watts\n', n_degree, predicted_metabolism_poly);
+        
+        % Plot results
+        figure;
+        plot(mass_kg, metabolism_watts, 'o', 'DisplayName', 'Data');
+        hold on;
+        x_plot = linspace(min(mass_kg), max(mass_kg) + 50, 100); % Extend plot range
+        y_fitted_poly = polyval(coefficients_poly, x_plot);
+        plot(x_plot, y_fitted_poly, 'b-', 'DisplayName', ['Polynomial Fit (Degree ' num2str(n_degree) ')']);
+        plot(tiger_mass, predicted_metabolism_poly, 'rx', 'MarkerSize', 10, 'LineWidth', 2, 'DisplayName', 'Tiger Prediction');
+        xlabel('Mass (kg)');
+        ylabel('Metabolism (watts)');
+        title(['Metabolism vs. Mass (polyfit, Degree ' num2str(n_degree) ')']);
+        legend('show');
+        grid on;
+        hold off;
+        
+    catch ME
+        fprintf('  Error encountered for degree %d: %s\n', n_degree, ME.message);
+        fprintf('  This degree might lead to a badly conditioned polynomial or a non-unique solution.\n');
+    end
+    
+    warning('on', 'MATLAB:polyfit:RepeatedRoots');
+    warning('on', 'MATLAB:polyfit:PolyNotUnique');
+end
+
+fprintf('\n  Comments on bad conditioning and non-unique fits for polyfit:\n');
+fprintf('  - **Badly conditioned:** This often occurs when fitting high-degree polynomials to a relatively small number of data points.\n');
+fprintf('    The system of equations that `polyfit` solves becomes sensitive to small changes in the data or numerical precision issues.\n');
+fprintf('    MATLAB might issue a warning (e.g., "Polynomial is badly conditioned"). This typically happens when the degree of the polynomial\n');
+fprintf('    approaches or exceeds the number of data points minus one.\n');
+fprintf('  - **No longer unique:** A polynomial fit of degree `n` to `m` data points is unique if and only if `n < m`. \n');
+fprintf('    If `n = m-1` (e.g., fitting a 5th-degree polynomial to 6 points), the polynomial will pass exactly through all points,\n');
+fprintf('    resulting in a unique polynomial, but it might oscillate wildly between points (overfitting).\n');
+fprintf('    If `n >= m`, the problem is ill-posed; infinitely many polynomials of degree `n` can pass through `m` points.\n');
+fprintf('    In your case, you have 6 data points. A unique fit is guaranteed for degrees up to 5. Trying degree 6 or higher \n');
+fprintf('    will lead to non-unique solutions or errors, as you''re trying to fit a higher-degree polynomial than the number of points allows.\n');
+
+
+%% (d) Do the tasks cited in (c) by using the Curve fitting Apps.
+
+fprintf('\n--- 1. Linear Regression (Curve Fitting App for Polynomials) ---\n');
+fprintf('  To perform this:\n');
+fprintf('  1. Go to the "APPS" tab -> "Curve Fitting".\n');
+fprintf('  2. Select Data: X=mass_kg, Y=metabolism_watts.\n');
+fprintf('  3. For "Fit type", select "Polynomial".\n');
+fprintf('  4. Change the "Degree" dropdown to try different polynomial degrees (e.g., 1, 2, 3, 4, 5).\n');
+fprintf('     Observe the fit and the coefficients. The app will visually show overfitting for high degrees.\n');
+fprintf('     It will also indicate the "Adjusted R-square" and "RMSE" for goodness of fit.\n');
+fprintf('  5. You can evaluate predictions for 350kg in the "Evaluate" panel.\n');
+fprintf('  Comments: The app provides a visual and interactive way to observe the effects of polynomial degree,\n');
+fprintf('  including potential overfitting and how R-squared values change.\n');
+
+
+%%% 2. Nonlinear Regression
+
+% Data
+P_obs = [0.985; 1.108; 1.363; 1.631]; % Pressure (P)
+V_obs = [25000; 22200; 18000; 15000]; % Volume (V)
+
+% Constants
+R = 82.05; % mL.atm/gmol.K
+T = 303;   % K
+
+% 1. Use the function fminsearch to write a script file using the following function file to find A1 and A2.
+% Make sure PVSSR.m is in a subfolder named 'lab4' or on MATLAB's path.
+% Initial guess for parameters [A1, A2]
+a0_fminsearch = [0.1; 0.1]; % Initial guess for A1 and A2
+
+% Define the objective function for fminsearch.
+% fminsearch minimizes a scalar function. Our PVSSR returns the sum of squared residuals.
+% So, we pass the parameters 'a' directly to PVSSR, along with the fixed x and y data.
+objective_func = @(a) PVSSR(a, V_obs, P_obs);
+
+% Perform the optimization
+[params_found_fminsearch, fval_fminsearch] = fminsearch(objective_func, a0_fminsearch);
+
+fprintf('\n--- 2. Nonlinear Regression (fminsearch) ---\n');
+fprintf('  Found parameters [A1, A2] using fminsearch: [ %.4f, %.4f ]\n', params_found_fminsearch(1), params_found_fminsearch(2));
+fprintf('  Minimum sum of squared residuals (fval): %.4e\n', fval_fminsearch);
+
+% Reconstruct and plot the fitted curve
+V_plot = linspace(min(V_obs)*0.9, max(V_obs)*1.1, 100); % Extended volume range for plotting
+P_fitted_fminsearch = R * T * (1 + params_found_fminsearch(1)./V_plot + params_found_fminsearch(2)./V_plot.^2) ./ V_plot;
+
+figure;
+plot(V_obs, P_obs, 'o', 'DisplayName', 'Observed Data');
+hold on;
+plot(V_plot, P_fitted_fminsearch, 'r-', 'DisplayName', 'Fitted Model');
+xlabel('Volume (mL)');
+ylabel('Pressure (atm)');
+title('Nonlinear Regression: Pressure vs. Volume');
+legend('show');
+grid on;
+hold off;
+
+
+% Explain: what is the meaning of a, x, y and f in PVSSR?
+fprintf('\n  Explanation of variables in PVSSR(a, x, y):\n');
+fprintf('  - **a**: This is a vector of the parameters that `fminsearch` is trying to find.\n');
+fprintf('           In this problem, `a` is `[A1, A2]`, where `a(1)` corresponds to `A1` and `a(2)` to `A2`.\n');
+fprintf('  - **x**: This is the independent variable data, which corresponds to `Volume (V)` in your problem.\n');
+fprintf('  - **y**: This is the dependent variable data, which corresponds to `Pressure (P)` in your problem.\n');
+fprintf('  - **f**: This is the output of the function, representing the **sum of squared residuals**.\n');
+fprintf('           It''s calculated as `sum((y_observed - y_predicted_by_model)^2)`. `fminsearch` tries to minimize this `f` value.\n');
+
+
+% 2. Use the function fminbnd to find the minimum value of the pressure if the volume varies between 2000 mL to 4000 mL.
+
+% The problem asks for the minimum value of Pressure (P) *as a function of V*,
+% using the A1 and A2 values we just found.
+% P(V) = R*T*(1 + A1/V + A2/V^2) / V
+% Define this as a function of V, using the found A1 and A2
+P_func_of_V = @(V) R * T * (1 + params_found_fminsearch(1)./V + params_found_fminsearch(2)./V.^2) ./ V;
+
+% Bounds for Volume
+V_lower_bound = 2000; % mL
+V_upper_bound = 4000; % mL
+
+% Find the minimum pressure within the given volume range
+[min_V_fminbnd, min_P_fminbnd] = fminbnd(P_func_of_V, V_lower_bound, V_upper_bound);
+
+fprintf('\n  Minimum Pressure (P) using fminbnd within V=[2000, 4000] mL:\n');
+fprintf('  Minimum Pressure = %.4f atm (at Volume = %.4f mL)\n', min_P_fminbnd, min_V_fminbnd);
+
+% Plot the pressure curve and mark the minimum
+figure;
+V_plot_bnd = linspace(V_lower_bound, V_upper_bound, 200);
+P_plot_bnd = P_func_of_V(V_plot_bnd);
+plot(V_plot_bnd, P_plot_bnd, 'b-', 'DisplayName', 'Pressure vs. Volume');
+hold on;
+plot(min_V_fminbnd, min_P_fminbnd, 'ro', 'MarkerSize', 8, 'LineWidth', 2, 'DisplayName', 'Minimum Pressure');
+xlabel('Volume (mL)');
+ylabel('Pressure (atm)');
+title('Pressure vs. Volume with Minimum (fminbnd)');
+legend('show');
+grid on;
+hold off;
+
+
+%%% 3. Descriptive Statistics in MATLAB
+
+%% a. Statistical functions in MATLAB
+
+% Student heights data (in cm)
+student_heights = [
+    176, 183, 152, 164, 159, 154;
+    167, 168, 156, 155, 153, 190;
+    182, 178, 159, 161, 166, 172;
+    174, 175, 168, 185, 184, 159;
+    164, 162, 164, 175, 170, 169
+];
+
+% Reshape to a single vector for easier processing
+s = student_heights(:); % Column vector of all heights
+
+fprintf('\n--- 3a. Descriptive Statistics (Student Heights) ---\n');
+
+% (a) the average height of the 30 students
+avg_height = mean(s);
+fprintf('  (a) Average height: %.2f cm\n', avg_height);
+
+% (b) the variance
+variance_height = var(s);
+fprintf('  (b) Variance: %.2f\n', variance_height);
+
+% (c) the standard deviation of the distribution
+std_dev_height = std(s);
+fprintf('  (c) Standard deviation: %.2f\n', std_dev_height);
+
+% (d) the lowest
+min_height = min(s);
+fprintf('  (d) Lowest height: %.0f cm\n', min_height);
+
+% (e) the highest heights of the student in this class
+max_height = max(s);
+fprintf('  (e) Highest height: %.0f cm\n', max_height);
+
+% (f) the mode
+mode_height = mode(s);
+fprintf('  (f) Mode height: %.0f cm\n', mode_height);
+
+% (g) the median height
+median_height = median(s);
+fprintf('  (g) Median height: %.0f cm\n', median_height);
+
+
+% Draw a histogram of the height of the students by dividing the students into 5 groups.
+num_bins = 5;
+figure;
+histogram(s, num_bins); % histogram is the modern function, hist is older
+% To get n and x for compatibility with the lab's hist description:
+[n_hist, x_hist] = hist(s, num_bins);
+fprintf('\n  Histogram Bin Counts (n): ');
+disp(n_hist'); % Transpose for column output
+fprintf('  Histogram Bin Midpoints (x): ');
+disp(x_hist'); % Transpose for column output
+
+xlabel('Height (cm)');
+ylabel('Number of Students');
+title(['Histogram of Student Heights (' num2str(num_bins) ' Groups)']);
+grid on;
+
+
+%% b. Uniform distribution and normal distribution in MATLAB
+
+fprintf('\n--- 3b. Uniform and Normal Distributions ---\n');
+
+% 100 - elements row vector A whose elements are randomly and uniformly distributed between 140 and 180
+len_A = 100;
+A = 140 + (180 - 140) * rand(1, len_A); % (b-a)*rand(size) + a
+
+fprintf('\n  --- Vector A (Uniform Distribution) ---\n');
+fprintf('  Size of A: %s\n', mat2str(size(A)));
+fprintf('  Average of A: %.4f\n', mean(A));
+fprintf('  Variance of A: %.4f\n', var(A));
+fprintf('  Standard Deviation of A: %.4f\n', std(A));
+fprintf('  Median of A: %.4f\n', median(A));
+fprintf('  Mode of A: %.4f (may not be very meaningful for continuous random data)\n', mode(A));
+% For continuous data, mode often returns the most frequent *exact* value, which is rare.
+% A histogram would be more informative for the mode of continuous data.
+
+
+% 200 - elements column vector B whose elements are randomly and normally distributed with mean 170 and standard deviation 30.
+len_B = 200;
+mean_B = 170;
+std_B = 30;
+B = mean_B + std_B * randn(len_B, 1); % mean + std * randn(size)
+
+fprintf('\n  --- Vector B (Normal Distribution) ---\n');
+fprintf('  Size of B: %s\n', mat2str(size(B)));
+fprintf('  Average of B: %.4f\n', mean(B));
+fprintf('  Variance of B: %.4f\n', var(B));
+fprintf('  Standard Deviation of B: %.4f\n', std(B));
+fprintf('  Median of B: %.4f\n', median(B));
+fprintf('  Mode of B: %.4f (may not be very meaningful for continuous random data)\n', mode(B));
+
+% Plot histograms for random data to visualize distributions
+figure;
+subplot(1,2,1);
+histogram(A);
+title('Histogram of Uniformly Distributed Data (A)');
+xlabel('Value');
+ylabel('Frequency');
+grid on;
+
+subplot(1,2,2);
+histogram(B);
+title('Histogram of Normally Distributed Data (B)');
+xlabel('Value');
+ylabel('Frequency');
+grid on;

linear_regression.m

@@ -0,0 +1,73 @@
+function [a0, a1, r2] = linear_regression(x, y)
+% linear_regression: Performs linear regression y = a0 + a1*x
+% x: independent variable (row vector)
+% y: dependent variable (row vector)
+% a0: y-intercept
+% a1: slope
+% r2: coefficient of determination
+
+if length(x) ~= length(y)
+    fprintf('The independent and the dependent variables should have same number of elements\n');
+    a0 = NaN; a1 = NaN; r2 = NaN; % Return NaNs if dimensions don't match
+    return; % Exit the function
+else
+    n = length(x); % The number of data
+end
+
+X = mean(x); Y = mean(y);
+
+% Ensure x and y are row vectors for dot products or use *'
+% The provided formula uses x*y' and x*x' which implies row vectors
+% if x and y are given as row vectors or column vectors and used with transpose.
+% Let's ensure they are row vectors for consistency with the formula.
+x = reshape(x, 1, []); % Ensure x is a row vector
+y = reshape(y, 1, []); % Ensure y is a row vector
+
+% The formula for a1 in the lab file: (x*y'/n - X*Y)/var(x);
+% This looks like sum(x.*y)/n - mean(x)*mean(y)
+% Let's clarify the formulas based on common linear regression.
+% a1 = (n*sum(x.*y) - sum(x)*sum(y)) / (n*sum(x.^2) - (sum(x)).^2);
+% Simplified, using deviations:
+
+Sxy = sum((x - X) .* (y - Y));
+Sxx = sum((x - X).^2);
+
+if Sxx == 0
+    % Handle case where all x values are the same (vertical line)
+    a1 = NaN;
+    a0 = NaN;
+    r2 = NaN;
+    fprintf('Cannot perform regression: all x values are the same.\n');
+    return;
+end
+
+a1 = Sxy / Sxx;
+a0 = Y - a1 * X;
+
+% R-squared calculation
+% r2= (n*x*y'-n*X*n*Y)/sqrt((n*x*x'-(n*X)^2)*(n*y*y'-(n*Y)^2));
+% This looks like the correlation coefficient 'r' squared.
+% The 'r' correlation coefficient: r = Sxy / sqrt(Sxx * Syy)
+% Where Syy = sum((y-Y).^2)
+
+Syy = sum((y - Y).^2);
+
+if Sxx == 0 || Syy == 0
+    r = NaN;
+else
+    r = Sxy / sqrt(Sxx * Syy);
+end
+r2 = r^2;
+
+% Create the best fit line for plotting
+xp = linspace(min(x), max(x), 2); % Only need 2 points for a straight line
+yp = a0 + a1 * xp;
+
+% This function also plots directly, which is useful for part (a)
+plot(x, y, 'o', xp, yp, 'b-');
+xlabel('Independent Variable (x)');
+ylabel('Dependent Variable (y)');
+title('Linear Regression Fit');
+grid on;
+
+end

newtmult.m

@@ -0,0 +1,66 @@
+% function [x,f,ea,iter]=newtmult(func,x0,es,maxit,varargin)
+% newtmult: Newton-Raphson root zeroes nonlinear systems
+% [x,f,ea,iter]=newtmult(func,x0,es,maxit,p1,p2,...):
+% uses the Newton-Raphson method to find the roots of
+% a system of nonlinear equations
+% input:
+% func = name of function that returns f and J
+% x0 = initial guess
+% es = desired percent relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by function
+% output:
+% x = vector of roots
+% f = vector of functions evaluated at roots
+% ea = approximate percent relative error (%)
+% iter = number of iterations
+% if nargin<2,error('at least 2 input arguments required'),end
+% if nargin<3||isempty(es),es=0.0001;end
+% if nargin<4||isempty(maxit),maxit=50;end
+% iter = 0;
+% x=x0;
+% while (1)
+% [J,f]=func(x,varargin{:});
+% dx=J\f;
+% x=x-dx;
+% a=100*max(abs(dx./x));
+% if iter>=maxit||ea<=es, break, end
+% end
+% end
+function [x, f, ea, iter] = newtmult(func, x0, es, maxit, varargin)
+% NEWTMULT: Newton-Raphson root zeroes nonlinear systems
+%   [x,f,ea,iter]=newtmult(func,x0,es,maxit,p1,p2,...)
+%   uses the Newton-Raphson method to find the roots of
+%   a system of nonlinear equations
+%
+%   input:
+%     func = name of function that returns f and J
+%     x0   = initial guess
+%     es   = desired percent relative error (default = 0.0001%)
+%     maxit= maximum allowable iterations (default = 50)
+%     p1,p2,... = additional parameters used by function
+%
+%   output:
+%     x    = vector of roots
+%     f    = vector of functions evaluated at roots
+%     ea   = approximate percent relative error (%)
+%     iter = number of iterations
+
+if nargin < 2, error('at least 2 input arguments required'), end
+if nargin < 3 || isempty(es), es = 0.0001; end
+if nargin < 4 || isempty(maxit), maxit = 50; end
+
+iter = 0;
+x = x0;
+
+while (1)
+    [J, f] = func(x, varargin{:}); % Key line: Get Jacobian and function values
+    dx = J \ f;                     % Solve for the change in x
+    x = x - dx;                     % Update x
+    iter = iter + 1;                % Increment iteration counter
+    ea = 100 * max(abs(dx ./ x));   % Calculate approximate relative error
+
+    if iter >= maxit || ea <= es, break; end % Check for convergence
+end
+
+end