restart;

with(LinearAlgebra): with(Statistics): with(plots):

local D:

# Coefficients for RSM model

beta0 := 1.5: beta1 := 0.8: beta2 := 0.5: beta3 := 0.7: beta4 := 0.3:

beta11 := 0.1: beta22 := 0.2: beta33 := 0.1: beta44 := 0.5e-1:

beta12 := 0.05: beta13 := 0.04: beta14 := 0.03: beta23 := 0.02:

beta24 := 0.03: beta34 := 0.02:

# RSM model for Nusselt number

Nu := beta0 + beta1*A + beta2*B + beta3*C + beta4*D +

      beta11*A^2 + beta22*B^2 + beta33*C^2 + beta44*D^2 +

      beta12*A*B + beta13*A*C + beta14*A*D +

      beta23*B*C + beta24*B*D + beta34*C*D:

# Central Composite Design (CCD)

DesignMatrix := Matrix([

[-1,-1,-1,-1],[1,-1,-1,-1],[-1,1,-1,-1],[1,1,-1,-1],

[-1,-1,1,-1],[1,-1,1,-1],[-1,1,1,-1],[1,1,1,-1],

[-1,-1,-1,1],[1,-1,-1,1],[-1,1,-1,1],[1,1,-1,1],

[-1,-1,1,1],[1,-1,1,1],[-1,1,1,1],[1,1,1,1],

[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]

]):

# Map coded levels to actual values

decode := proc(coded, low, high)

    return (high + low)/2 + coded*(high - low)/2;

end proc:

# Real value ranges

A_range := [0.5, 2.5]: B_range := [0.01, 0.05]:

C_range := [0.01, 0.05]: D_range := [0.0, 0.8]:

ResponseVector := Vector(20):

for i from 1 to 20 do

    A := decode(DesignMatrix[i,1], op(A_range));

    B := decode(DesignMatrix[i,2], op(B_range));

    C := decode(DesignMatrix[i,3], op(C_range));

    D := decode(DesignMatrix[i,4], op(D_range));

    ResponseVector[i] := evalf(eval(Nu), 8);

    printf("Run %d: [A = %.2f, B = %.2f, C = %.2f, D = %.2f] => Nu = %.8f\n", i, A, B, C, D, ResponseVector[i]);

end do:

pointplot([seq([i, ResponseVector[i]], i=1..20)],

          style=point, color=blue, symbol=circle,

          title="Nusselt Number vs. Run Index",

          labels=["Run Index", "Nusselt Number"]);

Run 1: [A = 0.50, B = 0.01, C = 0.01, D = 0.00] => Nu = 1.93748200
Run 2: [A = 2.50, B = 0.01, C = 0.01, D = 0.00] => Nu = 4.13928200
Run 3: [A = 0.50, B = 0.05, C = 0.01, D = 0.00] => Nu = 1.95897000
Run 4: [A = 2.50, B = 0.05, C = 0.01, D = 0.00] => Nu = 4.16477000
Run 5: [A = 0.50, B = 0.01, C = 0.05, D = 0.00] => Nu = 1.96653000
Run 6: [A = 2.50, B = 0.01, C = 0.05, D = 0.00] => Nu = 4.17153000
Run 7: [A = 0.50, B = 0.05, C = 0.05, D = 0.00] => Nu = 1.98805000
Run 8: [A = 2.50, B = 0.05, C = 0.05, D = 0.00] => Nu = 4.19705000
Run 9: [A = 0.50, B = 0.01, C = 0.01, D = 0.80] => Nu = 2.22188200
Run 10: [A = 2.50, B = 0.01, C = 0.01, D = 0.80] => Nu = 4.47168200
Run 11: [A = 0.50, B = 0.05, C = 0.01, D = 0.80] => Nu = 2.24433000
Run 12: [A = 2.50, B = 0.05, C = 0.01, D = 0.80] => Nu = 4.49813000
Run 13: [A = 0.50, B = 0.01, C = 0.05, D = 0.80] => Nu = 2.25157000
Run 14: [A = 2.50, B = 0.01, C = 0.05, D = 0.80] => Nu = 4.50457000
Run 15: [A = 0.50, B = 0.05, C = 0.05, D = 0.80] => Nu = 2.27405000
Run 16: [A = 2.50, B = 0.05, C = 0.05, D = 0.80] => Nu = 4.53105000
Run 17: [A = 1.50, B = 0.03, C = 0.03, D = 0.40] => Nu = 3.11193800
Run 18: [A = 1.50, B = 0.03, C = 0.03, D = 0.40] => Nu = 3.11193800
Run 19: [A = 1.50, B = 0.03, C = 0.03, D = 0.40] => Nu = 3.11193800
Run 20: [A = 1.50, B = 0.03, C = 0.03, D = 0.40] => Nu = 3.11193800