lda/nmu-m
2
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

added source file because of gitea webview not supporting .mlx

Browse Source 
added with save as
This commit is contained in:
lda committed 2025-06-14 03:58:12 +00:00
1 parent bb5767833b
commit 6e7b1062c1
5 files changed
+1549
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

lab5/lab5_t.m
+434
View File
@@ -0,0 +1,434 @@
%% Labwork 5: Linear Programming, Numerical Integration, Symbolic Math, Numerical Differentiation
%% 1. Linear Programming
%% Exercise 1: Gas Processing Plant Maximization Problem
% (a) Develop a linear programming formulation to maximize the profits for this operation.
% Let x1 = quantity of regular gas produced (in tons)
% Let x2 = quantity of premium gas produced (in tons)
% Objective Function (Maximize Profit):
% Max P = 150*x1 + 175*x2
% Constraints:
% 1. Raw Gas Constraint: 7*x1 + 11*x2 <= 77 (m^3/week)
% 2. Production Time Constraint: 10*x1 + 8*x2 <= 80 (hr/week)
% 3. Regular Gas Storage Constraint: x1 <= 9 (tons)
% 4. Premium Gas Storage Constraint: x2 <= 6 (tons)
% 5. Non-negativity Constraints: x1 >= 0, x2 >= 0
% (b) Minimization of the problem above using function: linprog in Matlab
% linprog solves minimization problems: min (f'*x) such that A*x <= b, Aeq*x = beq, lb <= x <= ub
% To convert our maximization problem (Max P = 150*x1 + 175*x2) to minimization:
% Min -P = -(150*x1 + 175*x2)
f = [-150; -175]; % Coefficients of the objective function for minimization
% Inequality constraints (A*x <= b)
A = [7 11; % Raw gas
10 8]; % Production time
b = [77; 80]; % Resource availability
% No equality constraints (Aeq, beq are empty or not used)
Aeq = [];
beq = [];
% Lower bounds (lb) and Upper bounds (ub) for x1 and x2
lb = [0; 0]; % Non-negativity constraints
ub = [9; 6]; % Storage constraints (x1 <= 9, x2 <= 6)
% Solve the linear programming problem
options = optimoptions('linprog', 'Display', 'off'); % Suppress iteration display
[x_opt, fval_opt, exitflag, output] = linprog(f, A, b, Aeq, beq, lb, ub, options);
fprintf('--- 1. Linear Programming (linprog) ---\n');
fprintf(' Optimal production of Regular Gas (x1): %.4f tons\n', x_opt(1));
fprintf(' Optimal production of Premium Gas (x2): %.4f tons\n', x_opt(2));
fprintf(' Maximum Profit: $%.2f\n', -fval_opt); % Negative of fval_opt because we minimized -P
% (c) Demonstrate the feasible solution space using the graphical solution of a linear programming problem.
% Hint: use linspace and meshgrid and plot
figure;
hold on;
xlabel('Regular Gas (x1)');
ylabel('Premium Gas (x2)');
title('Feasible Region and Optimal Solution for Gas Production');
grid on;
axis equal; % Ensures equal scaling on both axes
% Define range for plotting
x1_plot = linspace(0, 10, 200); % x1 from 0 to slightly above max possible
x2_plot = linspace(0, 10, 200); % x2 from 0 to slightly above max possible
% Plot constraints as lines (and shade feasible region)
% Constraint 1: 7*x1 + 11*x2 <= 77 => x2 <= (77 - 7*x1) / 11
plot(x1_plot, (77 - 7*x1_plot) / 11, 'DisplayName', 'Raw Gas Constraint (7x1 + 11x2 <= 77)', 'LineWidth', 1.5);
% Constraint 2: 10*x1 + 8*x2 <= 80 => x2 <= (80 - 10*x1) / 8
plot(x1_plot, (80 - 10*x1_plot) / 8, 'DisplayName', 'Production Time Constraint (10x1 + 8x2 <= 80)', 'LineWidth', 1.5);
% Constraint 3: x1 <= 9 (vertical line)
plot([9 9], [0 10], 'DisplayName', 'Storage Constraint (x1 <= 9)', 'LineWidth', 1.5);
% Constraint 4: x2 <= 6 (horizontal line)
plot([0 10], [6 6], 'DisplayName', 'Storage Constraint (x2 <= 6)', 'LineWidth', 1.5);
% Non-negativity constraints are implicit by starting x1_plot and x2_plot from 0
% Shade the feasible region (more complex for general cases, but for this 2D, we can draw a polygon)
% Find vertices of the feasible region
% (0,0)
% Intersection of x1=0 and x2=6 => (0,6)
% Intersection of x2=6 and 7x1 + 11x2 = 77 => 7x1 + 66 = 77 => 7x1 = 11 => x1 = 11/7 approx 1.57, (1.57, 6)
% Intersection of 7x1 + 11x2 = 77 and 10x1 + 8x2 = 80
% Multiply first by 8, second by 11:
% 56x1 + 88x2 = 616
% 110x1 + 88x2 = 880
% Subtracting first from second: 54x1 = 264 => x1 = 264/54 = 4.8889
% Then 10*(264/54) + 8x2 = 80 => 8x2 = 80 - 2640/54 = (4320-2640)/54 = 1680/54 => x2 = 1680/(54*8) = 1680/432 = 3.8889
% So, (4.8889, 3.8889)
% Intersection of 10x1 + 8x2 = 80 and x2=0 => x1 = 8
% Intersection of x1=9 and x2=0 => (9,0)
% Vertices to consider for feasible region:
% (0,0), (0, 6), (11/7, 6), (4.8889, 3.8889), (8, 0), (9, 0)
% We need to find the intersection of all binding constraints.
% The relevant vertices forming the feasible polygon (convex hull) are:
% (0,0)
% (0, 6) -- intersection of x1=0 and x2<=6
% (11/7, 6) ~= (1.5714, 6) -- intersection of 7x1 + 11x2 = 77 and x2=6
% (4.8889, 3.8889) -- intersection of 7x1 + 11x2 = 77 and 10x1 + 8x2 = 80
% (8, 0) -- intersection of 10x1 + 8x2 = 80 and x2=0
% (9, 0) -- intersection of x1=9 and x2=0
% It looks like the effective polygon vertices are (0,0), (0, 6), (11/7, 6), (4.8889, 3.8889), (8, 0)
feasible_x = [0, 0, 11/7, 264/54, 8];
feasible_y = [0, 6, 6, 1680/432, 0];
fill(feasible_x, feasible_y, [0.8 0.9 1], 'DisplayName', 'Feasible Region', 'FaceAlpha', 0.5); % Light blue
% Plot the optimal solution found by linprog
plot(x_opt(1), x_opt(2), 'ro', 'MarkerSize', 10, 'LineWidth', 2, 'MarkerFaceColor', 'r', 'DisplayName', 'Optimal Solution');
text(x_opt(1) + 0.2, x_opt(2) + 0.2, sprintf('(%.2f, %.2f)', x_opt(1), x_opt(2)), 'Color', 'r');
xlim([0 9.5]); % Adjust axis limits to focus on feasible region
ylim([0 7]);
legend('Location', 'best');
hold off;
%% 2. Numerical Integration
% Function definition
% fx = 0.2 + 25x - 200x2 + 675x3 - 900x4 + 400x5
f_integrand = @(x) 0.2 + 25*x - 200*x.^2 + 675*x.^3 - 900*x.^4 + 400*x.^5;
% Integration limits
a_int = 0;
b_int = 0.8;
% For trapz and cumtrapz, we need discrete points
num_points = 100; % Number of points for discretization
x_discrete = linspace(a_int, b_int, num_points);
y_discrete = f_integrand(x_discrete);
fprintf('\n--- 2. Numerical Integration ---\n');
%% (a) Use Matlab trapz function to estimate the integral above from 0 to 0.8
integral_trapz = trapz(x_discrete, y_discrete);
fprintf(' (a) Integral using trapz: %.8f\n', integral_trapz);
%% (b) Use the trapezoidal function in Matlab, cumtrapz. Compare the difference between two functions
% cumtrapz calculates the cumulative integral
cumulative_integral_cumtrapz = cumtrapz(x_discrete, y_discrete);
% The last value of cumtrapz should be equal to trapz for the whole interval
integral_cumtrapz_full = cumulative_integral_cumtrapz(end);
fprintf(' (b) Full integral using cumtrapz (last value): %.8f\n', integral_cumtrapz_full);
fprintf(' Difference between trapz and cumtrapz (full integral): %.8e\n', abs(integral_trapz - integral_cumtrapz_full));
fprintf(' Comment: trapz computes the total integral over the given range, while cumtrapz returns\n');
fprintf(' a vector where each element is the integral from the start point to that x-value.\n');
fprintf(' The last value of cumtrapz should conceptually match the result of trapz.\n');
%% (c) Use the simpson function in Matlab, quad
% The 'quad' function is largely superseded by 'integral' in modern MATLAB.
% 'integral' is more robust and generally preferred.
integral_quad_simpson = integral(f_integrand, a_int, b_int);
fprintf(' (c) Integral using integral (quad replacement): %.8f\n', integral_quad_simpson);
%% (d) Use the quadrature function quadadapt (custom function)
% Define the custom quadadapt and quadstep functions as local functions
% or ensure they are in separate .m files on your path.
% Function: q = quadadapt(f,a,b,tol,varargin)
% Recursive subfunction: q = quadstep(f,a,b,tol,fa,fc,fb,varargin)
% Custom quadadapt function (embed here for MLX, or ensure separate file)
% If this were a separate .m file, you'd only need the definition once.
% For MLX, define them inline.
% Make sure 'f' is a function handle.
integral_quadadapt = quadadapt(f_integrand, a_int, b_int);
fprintf(' (d) Integral using custom quadadapt: %.8f\n', integral_quadadapt);
% Embedded Custom Functions (for the .mlx file)
function q = quadadapt(f,a,b,tol,varargin)
% Evaluates definite integral of f(x) from a to b
if nargin < 4 || isempty(tol),tol = 1.e-6;end
c = (a + b)/2;
fa = feval(f,a,varargin{:});
fc = feval(f,c,varargin{:});
fb = feval(f,b,varargin{:});
q = quadstep(f, a, b, tol, fa, fc, fb, varargin{:});
end
function q = quadstep(f,a,b,tol,fa,fc,fb,varargin)
% Recursive subfunction used by quadadapt.
h = b - a; c = (a + b)/2;
fd = feval(f,(a+c)/2,varargin{:});
fe = feval(f,(c+b)/2,varargin{:});
q1 = h/6 * (fa + 4*fc + fb); % Simpson's Rule for the whole interval [a,b]
q2 = h/12 * (fa + 4*fd + 2*fc + 4*fe + fb); % Simpson's Rule for two sub-intervals [a,c] and [c,b]
if abs(q2 - q1) <= tol
q = q2 + (q2 - q1)/15; % Adaptive Simpson's correction
else
qa = quadstep(f, a, c, tol, fa, fd, fc, varargin{:});
qb = quadstep(f, c, b, tol, fc, fe, fb, varargin{:});
q = qa + qb;
end
end
%% 3. Symbolic Integration and Differentiation
fprintf('\n--- 3. Symbolic Integration and Differentiation ---\n');
%% a. Define the following function using syms: fx = x2 ex - 5x3
% Compute the integral, and the first and second derivatives of the above function
% symbolically
syms x_sym
f_sym_a = x_sym^2 * exp(x_sym) - 5*x_sym^3;
% Compute the integral
integral_f_sym_a = int(f_sym_a, x_sym);
% Compute the first derivative
first_deriv_f_sym_a = diff(f_sym_a, x_sym);
% Compute the second derivative
second_deriv_f_sym_a = diff(f_sym_a, x_sym, 2);
fprintf(' (a) Original function f(x): %s\n', char(f_sym_a));
fprintf(' Symbolic Integral of f(x): %s\n', char(integral_f_sym_a));
fprintf(' Symbolic First Derivative of f(x): %s\n', char(first_deriv_f_sym_a));
fprintf(' Symbolic Second Derivative of f(x): %s\n', char(second_deriv_f_sym_a));
%% b. Define the following function using syms: fxyz = x2 ey - 5z2
% Compute the integral with respect to x and second derivative with respect
% to z
syms x_sym y_sym z_sym
f_sym_b = x_sym^2 * exp(y_sym) - 5*z_sym^2;
% Compute the integral with respect to x
integral_f_sym_b_x = int(f_sym_b, x_sym);
% Compute the second derivative with respect to z
second_deriv_f_sym_b_z = diff(f_sym_b, z_sym, 2);
fprintf('\n (b) Original function f(x,y,z): %s\n', char(f_sym_b));
fprintf(' Symbolic Integral of f(x,y,z) with respect to x: %s\n', char(integral_f_sym_b_x));
fprintf(' Symbolic Second Derivative of f(x,y,z) with respect to z: %s\n', char(second_deriv_f_sym_b_z));
%% 4. Numerical Differentiation
% Re-use the function f(x) from Exercise 2 for numerical differentiation
% fx = 0.2 + 25x - 200x2 + 675x3 - 900x4 + 400x5
f_num_diff = @(x) 0.2 + 25*x - 200*x.^2 + 675*x.^3 - 900*x.^4 + 400*x.^5;
% Define a range for x to calculate derivatives
x_range_diff = 0:0.01:0.8; % A finer step for numerical differentiation
y_range_diff = f_num_diff(x_range_diff);
fprintf('\n--- 4. Numerical Differentiation ---\n');
%% a. Use Matlab diff function to differentiate the function f(x) from Ex2
% When used on arrays, diff calculates discrete differences. dy/dx approx diff(y)./diff(x)
% First derivative
dydx_numerical = diff(y_range_diff) ./ diff(x_range_diff);
% The resulting derivative array is one element shorter than the original array.
% The x-values for the derivative are typically taken as the midpoints of the original x-intervals.
x_deriv_plot = (x_range_diff(1:end-1) + x_range_diff(2:end))/2;
figure;
plot(x_deriv_plot, dydx_numerical, 'b-', 'DisplayName', 'Numerical Derivative (diff)');
hold on;
xlabel('x');
ylabel('f''(x)');
title('Numerical Differentiation using diff()');
grid on;
%% b. Comparison of the exact derivative with numerical estimates
syms x_sym_diff
f_sym_diff_exact = 0.2 + 25*x_sym_diff - 200*x_sym_diff^2 + 675*x_sym_diff^3 - 900*x_sym_diff^4 + 400*x_sym_diff^5;
exact_deriv_f_sym = diff(f_sym_diff_exact, x_sym_diff);
% Convert symbolic derivative to a numeric function handle
exact_deriv_func = matlabFunction(exact_deriv_f_sym);
% Calculate exact derivative values over the same x range
exact_deriv_values = exact_deriv_func(x_deriv_plot); % Use x_deriv_plot for comparison
plot(x_deriv_plot, exact_deriv_values, 'r--', 'DisplayName', 'Exact Derivative');
legend('show');
hold off;
% Calculate average absolute difference
avg_abs_diff = mean(abs(dydx_numerical - exact_deriv_values));
fprintf(' (b) Average absolute difference between numerical (diff) and exact derivative: %.4e\n', avg_abs_diff);
%% c. Use forward function
% y = vector containing function values x = vector containing increments p =
% point number where to calculate d = derivative of y at point p Custom forward
% difference function (embed here or ensure separate file) function d = forward(y,x,p)
% d = (y(p+1)-y(p))/(x(p+1)-x(p)); end
% Test a point, e.g., at x=0.4 (which is at index 41 if step is 0.01 from 0)
test_point_index = find(x_range_diff == 0.4);
if isempty(test_point_index) || test_point_index == length(x_range_diff)
test_point_index = floor(length(x_range_diff)/2); % Fallback to middle if 0.4 not exact or is last point
fprintf(' Warning: Exact x=0.4 not found or is last point for forward diff. Using index %d.\n', test_point_index);
end
forward_deriv_at_point = forward(y_range_diff, x_range_diff, test_point_index);
exact_deriv_at_point = exact_deriv_func(x_range_diff(test_point_index));
fprintf('\n --- Numerical Differentiation (Custom Methods) ---\n');
fprintf(' (c) Forward difference at x=%.2f (index %d): %.8f (Exact: %.8f)\n', ...
x_range_diff(test_point_index), test_point_index, forward_deriv_at_point, exact_deriv_at_point);
%% d. Use backward function
% function d = backward(y,x,p) d = (y(p)-y(p-1))/(x(p)-x(p-1)); end
% Test a point, e.g., at x=0.4 (index 41)
% Ensure p is not the first element
if test_point_index == 1
test_point_index_bwd = 2; % Use second element for backward diff
else
test_point_index_bwd = test_point_index;
end
backward_deriv_at_point = backward(y_range_diff, x_range_diff, test_point_index_bwd);
exact_deriv_at_point_bwd = exact_deriv_func(x_range_diff(test_point_index_bwd));
fprintf(' (d) Backward difference at x=%.2f (index %d): %.8f (Exact: %.8f)\n', ...
x_range_diff(test_point_index_bwd), test_point_index_bwd, backward_deriv_at_point, exact_deriv_at_point_bwd);
%% e. Use central function
% function d = central(y,x,p) d = (y(p+1)-y(p-1))/(x(p+1)-x(p-1)); end
% Test a point, e.g., at x=0.4 (index 41)
% Ensure p is not the first or last element
if test_point_index == 1
test_point_index_central = 2; % Use second element for central diff
elseif test_point_index == length(x_range_diff)
test_point_index_central = length(x_range_diff) - 1; % Use second to last element
else
test_point_index_central = test_point_index;
end
central_deriv_at_point = central(y_range_diff, x_range_diff, test_point_index_central);
exact_deriv_at_point_central = exact_deriv_func(x_range_diff(test_point_index_central));
fprintf(' (e) Central difference at x=%.2f (index %d): %.8f (Exact: %.8f)\n', ...
x_range_diff(test_point_index_central), test_point_index_central, central_deriv_at_point, exact_deriv_at_point_central);
fprintf('\n Comment: Central difference formula generally provides a more accurate approximation\n');
fprintf(' of the derivative than forward or backward differences for a given step size.\n');
% --- Embedded Custom Functions (for the .mlx file) ---
% These functions are defined here to be self-contained within the .mlx
% For production code, it's usually better to have them in separate .m files.
% Custom linear_regression function (from Lab 4B, Ex 1a)
function [a0, a1, r2] = linear_regression(x, y)
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;
else
n = length(x);
end
X = mean(x); Y = mean(y);
x = reshape(x, 1, []);
y = reshape(y, 1, []);
Sxy = sum((x - X) .* (y - Y));
Sxx = sum((x - X).^2);
if Sxx == 0
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;
Syy = sum((y - Y).^2);
if Sxx == 0 || Syy == 0
r = NaN;
else
r = Sxy / sqrt(Sxx * Syy);
end
r2 = r^2;
xp = linspace(min(x), max(x), 2);
yp = a0 + a1 * xp;
figure; % Create new figure for this plot
plot(x, y, 'o', xp, yp, 'b-');
xlabel('Independent Variable (x)');
ylabel('Dependent Variable (y)');
title('Linear Regression Fit (from linear_regression.m)');
grid on;
end
% Custom PVSSR function (from Lab 4B, Ex 2)
% NOTE: This function is expected to be in a subfolder like './lab4/PVSSR.m'
% For this MLX, I'll include its code definition as if it were there.
% If you actually have it in a separate file, ensure MATLAB's path includes it.
% function f = PVSSR(a, x, y)
% R = 82.05;
% T = 303;
%
% if length(x) ~= length(y)
% fprintf('The two data must have the same length\n');
% f = NaN;
% return;
% end
%
% n = length(x);
%
% yp_model = zeros(1, n);
% for j = 1:n
% yp_model(j) = R * T * (1 + a(1)/x(j) + a(2)/x(j)^2) / x(j);
% end
%
% y = reshape(y, 1, []);
% yp_model = reshape(yp_model, 1, []);
%
% f = sum((y - yp_model).^2);
% end
% this is now at ./ as it is a common function
% Custom forward difference function (from Lab 5, Ex 4c)
function d = forward(y,x,p)
d = (y(p+1)-y(p))/(x(p+1)-x(p));
end
% Custom backward difference function (from Lab 5, Ex 4d)
function d = backward(y,x,p)
d = (y(p)-y(p-1))/(x(p)-x(p-1));
end
% Custom central difference function (from Lab 5, Ex 4e)
function d = central(y,x,p)
d = (y(p+1)-y(p-1))/(x(p+1)-x(p-1));
end