QUESTION 5

import numpy as np

import [Link] as plt

from scipy import stats

# Given Data

t = [Link]([0, 2, 4, 6, 8])

Ca = [Link]([1, 0.82, 0.67, 0.55, 0.45])

# Linearize the data: ln(Ca) = ln(Ca0) - k*t

ln_Ca = [Link](Ca)

# Perform linear least squares regression

slope, intercept, r_value, p_value, std_err = [Link](t, ln_Ca)

# Extract rate constant (k) and R-squared

k = -slope

R_squared = r_value**2

Ca0_calc = [Link](intercept)

# Verification: Calculate predicted values using the obtained k and Ca0_calc

Ca_pred = Ca0_calc * [Link](-k * t)

# Calculate percentage errors

errors = [Link]((Ca - Ca_pred) / Ca) * 100

print(f"Rate constant (k) = {k:.5f} min^-1")

print(f"Calculated Ca0 = {Ca0_calc:.5f} mol/L")

print(f"Coefficient of determination (R^2) = {R_squared:.5f}\n")

Rate constant (k) = 0.09982 min^-1

Calculated Ca0 = 1.00023 mol/L

Coefficient of determination (R^2) = 0.99999

import numpy as np

def solve_cstr():

v=1

V = 100

k1 = 1

k2 = 1

CA0 = 1

CB0 = 2

CC0 = 0

CD0 = 0

def F(x):

x1, x2, x3, x4 = x

f1 = v * (CA0 - x1) - V * k1 * x1 * x2

f2 = v * (CB0 - x2) - V * (k1 * x1 * x2 + k2 * x3 * x2)

f3 = v * (CC0 - x3) + V * (k1 * x1 * x2 - k2 * x3 * x2)

f4 = v * (CD0 - x4) + V * k2 * x3 * x2

return [Link]([f1, f2, f3, f4])

def Jacobian(x):

x1, x2, x3, x4 = x

J = [Link]((4, 4))

# f1 derivatives

J[0, 0] = -v - V * k1 * x2

J[0, 1] = -V * k1 * x1

J[0, 2] = 0

J[0, 3] = 0
 # f2 derivatives

J[1, 0] = -V * k1 * x2

J[1, 1] = -v - V * k1 * x1 - V * k2 * x3

J[1, 2] = -V * k2 * x2

J[1, 3] = 0

# f3 derivatives

J[2, 0] = V * k1 * x2

J[2, 1] = V * k1 * x1 - V * k2 * x3

J[2, 2] = -v - V * k2 * x2

J[2, 3] = 0

# f4 derivatives

J[3, 0] = 0

J[3, 1] = V * k2 * x3

J[3, 2] = V * k2 * x2

J[3, 3] = -v

return J

x = [Link]([0.2, 0.2, 0.2, 0.5])

print(f"Initial Guess: x = {x}")

for i in range(1, 4):

f_val = F(x)

j_val = Jacobian(x)

delta = [Link](j_val, -f_val)

x = x + delta

print(f"Iteration {i}: x = {x}")

print(f"F(x) at iter {i}: {F(x)}")

RESULTS

Initial Guess: x = [0.2 0.2 0.2 0.5]

Iteration 1: x = [0.06246282 0.18441404 0.0594884 0.87804878]

F(x) at iter 1: [-0.21436491 -0.43336574 -0.00463592 0.21900083]

Iteration 2: x = [0.05688275 0.16746306 0.05369756 0.88941969]

F(x) at iter 2: [-0.00945876 -0.0192748 -0.00035729 0.00981604]

Iteration 3: x = [0.05661427 0.16663781 0.05340926 0.88997646]

F(x) at iter 3: [-2.21562215e-05 -4.59477883e-05 -1.63534542e-06 2.37915669e-05]

import numpy as np

T = [Link]([300, 400, 500, 600])

Cp = [Link]([45.2, 52.1, 61.5, 73.4])

# Fit to 2nd order polynomial: Cp = a0 + a1*T + a2*T^2

# [Link] returns coefficients from highest degree to lowest: [a2, a1, a0]

coeffs = [Link](T, Cp, 2)

a2, a1, a0 = coeffs

print(f"a0 = {a0:.5f}")

print(f"a1 = {a1:.5f}")

print(f"a2 = {a2:.5f}")

RESULTS

a0 = 39.50000

a1 = -0.01850

a2 = 0.00013

QUESTION 1

import numpy as np

# 1. Define the coefficient matrix A and constant vector b

A = [Link]([

[5.0, 2.0, 3.0],

[3.0, 5.0, 2.0],

[2.0, 3.0, 5.0]

])

b = [Link]([350.0, 350.0, 300.0])

def gauss_elimination(A, b):

n = len(b)

Aug = [Link]((A, [Link](-1, 1))) # Augmented matrix

# Forward Elimination

for i in range(n):

for j in range(i + 1, n):

factor = Aug[j, i] / Aug[i, i]

Aug[j] = Aug[j] - factor * Aug[i]

# Backward Substitution

x = [Link](n)

for i in range(n - 1, -1, -1):

x[i] = (Aug[i, -1] - [Link](Aug[i, i+1:n] * x[i+1:n])) / Aug[i, i]

return x, Aug

# 3. Execute and Print Results

flow_rates, final_matrix = gauss_elimination(A, b)

print("--- Final Upper Triangular Augmented Matrix ---")

print([Link](final_matrix, 4))

print("\n--- Required Feed Flow Rates (kmol/h) ---")

print(f"Feed Stream 1 (x1): {flow_rates[0]:.2f}")

print(f"Feed Stream 2 (x2): {flow_rates[1]:.2f}")

print(f"Feed Stream 3 (x3): {flow_rates[2]:.2f}")

RESULTS

Feed Stream 1 (x1): 42.86

Feed Stream 2 (x2): 35.71

QUESTION 3

import numpy as np

# 1. Define the reordered, diagonally dominant system

A = [Link]([

[5.0, -1.0, 0.0],

[-1.0, 5.0, -1.0],

[ 0.0, -1.0, 4.0]

])

b = [Link]([9.0, 7.0, 5.0])

# 2. Set Parameters

omega = 1.2

iterations = 3

x = [Link](3) # Initial guess [0, 0, 0]

# 3. SOR Algorithm

print("--- Successive Over-Relaxation (SOR) Iterations ---")

print(f"Iteration 0: x1 = {x[0]:.4f}, x2 = {x[1]:.4f}, x3 = {x[2]:.4f}")

for k in range(1, iterations + 1):

x_old = [Link]()

for i in range(len(b)):

# Calculate the sum of A_ij * x_j (excluding the diagonal element)

sigma = sum(A[i][j] * x[j] for j in range(len(b)) if j != i)

# Gauss-Seidel estimate

x_gs = (b[i] - sigma) / A[i][i]

x[i] = (1 - omega) * x_old[i] + omega * x_gs

print(f"Iteration {k}: x1 = {x[0]:.4f}, x2 = {x[1]:.4f}, x3 = {x[2]:.4f}")

print("\n--- Final Estimated Compositions ---")

print(f"Tray 1 (x1): {x[0]:.4f}")

print(f"Tray 2 (x2): {x[1]:.4f}")

print(f"Tray 3 (x3): {x[2]:.4f}")

RESULTS

--- Final Estimated Compositions ---

Tray 1 (x1): 2.2609

Tray 2 (x2): 2.1846

Tray 3 (x3): 1.8038

import numpy as np

# 1. Define the System Matrices

A = [Link]([

[4.2, -1.1, -0.8, 0.0],

[-1.2, 5.3, -1.5, -1.0],

[-0.5, -1.4, 4.8, -1.2],

[0.0, -0.9, -1.3, 4.5]

])

B = [Link]([150.0, 200.0, 180.0, 160.0])

# 2. LU Decomposition Function (Doolittle's Algorithm)

def lu_decomposition(A):

n = len(A)

L = [Link]((n, n))

U = [Link]((n, n))

for i in range(n):

for k in range(i, n):

sum_u = sum(L[i][j] * U[j][k] for j in range(i))

U[i][k] = A[i][k] - sum_u

for k in range(i, n):

if i == k:

L[i][i] = 1.0 # Diagonal as 1

else:

sum_l = sum(L[k][j] * U[j][i] for j in range(i))

L[k][i] = (A[k][i] - sum_l) / U[i][i]

return L, U
# 3. Forward and Backward Substitution

def forward_substitution(L, B):

n = len(B)

y = [Link](n)

for i in range(n):

y[i] = B[i] - sum(L[i][j] * y[j] for j in range(i))

return y

def backward_substitution(U, y):

n = len(y)

T = [Link](n)

for i in range(n - 1, -1, -1):

T[i] = (y[i] - sum(U[i][j] * T[j] for j in range(i + 1, n))) / U[i][i]

return T

# 4. Execute and Print Results

L, U = lu_decomposition(A)

y = forward_substitution(L, B)

T = backward_substitution(U, y)

np.set_printoptions(precision=4, suppress=True)

print("--- Lower Matrix (L) ---\n", L)

print("\n--- Upper Matrix (U) ---\n", U)

print("\n--- Intermediate Vector (y) ---\n", y)

print("\n--- Final Temperatures (°C) ---")

for i, temp in enumerate(T, 1):

print(f"T{i} = {temp:.2f}")

results

--- Lower Matrix (L) ---

[[ 1. 0. 0. 0. ]
[-0.2857 1. 0. 0. ]

[-0.119 -0.3071 1. 0. ]

[ 0. -0.1805 -0.3862 1. ]]

--- Upper Matrix (U) ---

[[ 4.2 -1.1 -0.8 0. ]

[ 0. 4.9857 -1.7286 -1. ]

[ 0. 0. 4.174 -1.5071]

[ 0. 0. 0. 3.7374]]

--- Intermediate Vector (y) ---

[150. 242.8571 272.4308 309.0553]

--- Final Temperatures (°C) ---

T1 = 79.57

T2 = 98.28

T3 = 95.13

T4 = 82.69

Question 6

import numpy as np

# 1. Input the experimental data

t = [Link]([0, 2, 4, 6, 8])

Ca = [Link]([1.00, 0.82, 0.67, 0.55, 0.45])

# 2. Linearize: ln(Ca)

ln_Ca = [Link](Ca)

# 3. Perform linear least squares regression manually

# Calculate slope (m) and intercept (c)

sum_x = [Link](t)

sum_y = [Link](ln_Ca)

sum_xy = [Link](t * ln_Ca)

sum_x2 = [Link](t**2)

slope = (n * sum_xy - sum_x * sum_y) / (n * sum_x2 - sum_x**2)

intercept = (sum_y - slope * sum_x) / n

# 4. Calculate rate constant and R^2

k = -slope

# Predicted values

ln_Ca_pred = slope * t + intercept

# R-squared calculation

ss_total = [Link]((ln_Ca - [Link](ln_Ca))**2)

ss_res = [Link]((ln_Ca - ln_Ca_pred)**2)

R_squared = 1 - (ss_res / ss_total)

# 5. Display results

print("--- Linear Least Squares Regression Results ---")

print(f"Slope (m) = {slope:.4f}")

print(f"y-intercept (c) = {intercept:.4f}")

print("-" * 45)

print(f"Rate constant (k) = {k:.4f} min^-1")

print(f"R-squared (R^2) = {R_squared:.4f}")

RESULTS
Slope (m) = -0.0998

y-intercept (c) = 0.0002

---------------------------------------------

Rate constant (k) = 0.0998 min^-1

R-squared (R^2) = 0.9999