restart;

kernelopts(version);

A boundary value problem for an ODE

Let's solve the following boundary value problem:

Maple's dsolve() fails to find a solution:

Solution obtained by hand

The solution may be obtained by hand, as follows.

Integrate the DE once

.

where  is the integration constant.  Therefore

.

To determine the value of the constant , integrate the above over

and apply the boundary conditions:

whence
.

Having thus determined , we integrate the expression for  obtained above, and arrive at the solution
.

Remark: The solution obtained here is valid for any  as long as

the integrals encountered above make sense.  Note that there is
no requirement on the differentiability or even continuity of .

Technically,  can be any function such that  is integrable.

Wrong solution returned by dsolve()

Let us consider a special choice of the coefficient :

According the the calculations above, we have:

and therefore the correct solution is

Here is what the solution looks like:

Let's attempt to solve the same boundary value problem with Maple's dsolve().

That's very different from the correct solution obtained above.   To see the problem with it, let's recall that according to the DE, the expression should evaluate to   for some constant ,  but what we get through bad_sol is nothing like ;  in fact, it's not even continuous:

As solution that you  provided is a weak solution of the boundary value problem, (lets sayBVP(0)) one might think to it as the limit solution of
some regularized ("viscous") boundary value problem BVP(λ), λ ≥ 0, when  (from the right).
With this in mind on may think to  as some "shock" ad a way to damp it by adding a "viscosity" can be modeled by replacing  by :

The limit case  (from the right)
Then

Case x < 0

Case x > 0

restart

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

local A, B, C, D;

Digits := 10:

# Step 1: Define the variables

vars := [A, B, C, D];

# Step 2: Generate Central Composite Design (CCD) coded matrix

coded := 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]

]):

# Step 3: Define uncoded variable ranges

rangeA := [0.1, 0.5]:

rangeB := [0.1, 0.5]:

rangeC := [0.1, 0.5]:

rangeD := [0.1, 0.5]:

decode := (val, r) -> r[1] + (val + 1)*(r[2] - r[1])/2:

# Step 4: Create actual values matrix

actual := Matrix(20, 4):

for i from 1 to 20 do

    for j from 1 to 4 do

        r := eval(cat("range", vars[j])):

        actual[i, j] := evalf(decode(coded[i, j], r));

    end do;

end do:


# Verify "actual":

actual[1..2];

# I propose you this construction (it gives the correct result and avoid loops)

Low  := DiagonalMatrix(map2(op, 1, [rangeA, rangeB, rangeC, rangeD])):
High := DiagonalMatrix(map2(op, 2, [rangeA, rangeB, rangeC, rangeD])):

actual := coded . (High - Low) /~ 2 + Matrix(20, 4, 1) . (High + Low) /~ 2;

# Step 5: Define model coefficients (for simulation purposes)

beta := [1.5, 0.8, 0.6, 0.5, 0.3,  # Linear terms

         0.1, 0.2, 0.15, 0.18,  # Square terms

         0.05, 0.04, 0.03, 0.02, 0.01, 0.025]; # Interaction terms
vars

# Step 6: Simulate Response (Nusselt number)
Nu := Vector(20):

for i from 1 to 20 do
    # The next line, if executed changes changes "var".
    # indeed A is set to O.3 and thus vars = [A, B, C, D] becomes [0.3, B, C, D] and so on

    # A := actual[i,1]; B := actual[i,2]; C := actual[i,3]; D := actual[i,4];

    Nu[i] := evalf(
               eval(

                 beta[1] + beta[2]*A + beta[3]*B + beta[4]*C + beta[5]*D +

                 beta[6]*A^2 + beta[7]*B^2 + beta[8]*C^2 + beta[9]*D^2 +

                 beta[10]*A*B + beta[11]*A*C + beta[12]*A*D +

                 beta[13]*B*C + beta[14]*B*D + beta[15]*C*D
                 ,
                 vars =~ convert(actual[i], list)
               )

             ):

end do:

# Verify that "Nu" is a numeric vector

Nu[1..5];

# Verify that "vars" has not been changed

vars

# My Step 7a: Build the coded design matrix

basis           := [ 1, A, B, C, D, A^2, B^2, C^2, D^2, A*B, A*C, A*D, B*C, B*D, C*D]:
StructuralModel := add(cat(`alpha__`, b) * b, b in basis);

N := numelems(coded[..,1]);

CodedDesign := `<|>`( Vector(N, 1), coded ):

for i from 6 to numelems(basis) do
  opb := op(basis[i]):
  if opb[2]::numeric then
    col := ListTools:-Search(opb[1], basis) - 1;
    CodedDesign := `<|>`( CodedDesign, coded[.., col]^~opb[2] );
  else
    col := ListTools:-Search(opb[1], basis) - 1, ListTools:-Search(opb[2], basis) - 1;
    CodedDesign := `<|>`( CodedDesign, coded[.., col[1]]*~coded[.., col[2]] ):
  end if;
end do:

CodedDesign;
 

# My Step 7b: Compute the rank of "CodedDesign"

Rank(CodedDesign);

# You observe that "CodedDesign" is not of full rank (15) which means that it's not an orthogonal
# design matrix and thus some effects are aliased with some others.
#

# My Step 7c: Aliased effets:

for i from 1 to numelems(basis)-1 do
  for j from i+1 to numelems(basis) do
    if [min, max](CodedDesign[.., i] - CodedDesign[.., j]) = [0, 0] then
      print(cat("Aliasing between effects ", basis[i], " and ", basis[j]))
    end if:
  end do:
end do:

# So all the quadratic terms are mutualy aliased and cannot be identified independently
# Moer precisely it is impossible to say is some observed quadratic effect is A^2, or
# B^2, or C^2, or D^2

# My Step 7b: Build the design matrix with quadratic terms removed

reduced_basis   := [ 1, A, B, C, D, A*B, A*C, A*D, B*C, B*D, C*D]:
StructuralModel := add(cat(`alpha__`, b) * b, b in reduced_basis);

N := numelems(actual[..,1]);

DesignMatrix := `<|>`( Vector(N, 1), actual );
for i from 6 to numelems(reduced_basis) do
  opb := op(reduced_basis[i]):
  col := ListTools:-Search(opb[1], reduced_basis) - 1, ListTools:-Search(opb[2], reduced_basis) - 1;
  DesignMatrix := `<|>`( DesignMatrix, actual[.., col[1]]*~actual[.., col[2]] ):
  print(i, reduced_basis[i], [opb], [col]):
end do:

DesignMatrix;

Rank(DesignMatrix); # should be numelems(reduced_basis), thus 11
 

# The direct way to get the coefficient of the model

model_coeffs := (DesignMatrix^+ . DesignMatrix)^(-1) . (DesignMatrix^+ . Nu);

model_fitted := < reduced_basis > . model_coeffs assuming real;
 

# compare to your a prior model

evalf[3](model_fitted);

model_prior := beta[1] + beta[2]*A + beta[3]*B + beta[4]*C + beta[5]*D +
               beta[6]*A^2 + beta[7]*B^2 + beta[8]*C^2 + beta[9]*D^2 +
               beta[10]*A*B + beta[11]*A*C + beta[12]*A*D +
               beta[13]*B*C + beta[14]*B*D + beta[15]*C*D;


cmf := coeffs(evalf[3](model_fitted), vars, 'fmon');
fmon;

cmp := coeffs(model_prior, vars, 'cmon');
cmon;

# Agreement between "beta" and "cmf"

interface(rtablesize=16):
`<,>`(
  `<|>`(["effect", "beta", "fit"]),
  `<|>`(
    `<,>`(fmon),
    `<,>`(cmp[1..numelems([fmon])]),
    `<,>`(cmf)
  )
);

# the discrepancy between columns 2 and 3 comes from the quadratic effects which are considered
# as a model error or as an "observation" error due to some additive "noise.

# Step 8: Using Fit
model, Alpha := Fit(
           StructuralModel,
           actual,
           Nu,
           vars,
           output=["leastsquaresfunction", "parametervalues"]
         )[]:

model;

print():

# Observe you get the same results as above).

evalf[3](Alpha)

# Finally

Predictor := unapply(model_fitted, vars);

display(
  pointplot([seq([i, Nu[i]], i = 1..20)], style=point, symbol=circle,

          color=red, title="Nusselt Number vs Run Index",

          labels=["Run", "Nusselt Number"]
          , legend="Observations")
  ,
  pointplot([seq([i, Predictor(convert(actual[i], list)[])], i = 1..20)], style=point, symbol=cross, color=blue
          , legend="Predictions")
);
  

# I don't understand what you are trying to achieve using ANOVA?
#
# Clarification is required

# Step 9: ANOVA table

anovaTable := ANOVA(model, significancelevel=0.05):

# Step 10: Display results with 8 decimal places

printf("RSM Model for Nu (Nusselt number):\n");

print(evalf[8](model));

printf("\nANOVA Table:\n");

print(anovaTable);

restart

approx_sol := 3 -8*x + 24*x^2 - 208/3*x^3 + 200*x^4 - 8656/15*x^5 + 24976/15*x^6 + 553339/315*x^7 + 15550/21*x^8 + 502784/2835*x^9 + 125536/4725*x^10 + 78362/31185*x^11

p0 := numapprox:-pade(approx_sol, x, [3,3]);

p1 := map(sort, p0, x, ascending)