Maximizing arbitrarily scaled step function

Question

I am attempting to maximize a function dependent on a single variable (x) that contains multiple parameters. The case of a simply loaded beam with discontinuities in the loading. However so far I have only been able to get an answer using findMaximum and real values for the parameters, is it possible to get an output in terms of some constant and P, L, and EI?

    Rab = 13 P/32;
    Rbb = 27 P/32;
    a = L/4;
    b = L/2;
    c = 3 L/4;
    q = P/L;
    load[x_] = 
    Rab DiracDelta[x] - q (HeavisideTheta[x - a] - HeavisideTheta[x - b]) - P DiracDelta[x - c] + Rbb DiracDelta[x - L];
    ode = EI v''''[x] == load[x];
    c = {v[0] == 0, v[L] == 0, v''[0] == 0, v''[L] == 0};
    out = DSolve[{ode, c}, v[x], x];
    v[x_] = v[x] /. out;
    s[x_] = v'[x];
    V[x_] = -EI v'''[x];
    M[x_] = EI v''[x];
    Maximize[V[x], x]

This code simply returns the input to the Maximize function, V[x].

Answer 1

Since differential equation has only fourth derivative on lhs, integrate four times and determine integration variables from initial conditions.

Rab = 13 P/32;
Rbb = 27 P/32;
a = L/4;
b = L/2;
c = 3 L/4;
q = P/L;
load[x_] = 
  Rab DiracDelta[x] - 
   q (HeavisideTheta[x - a] - HeavisideTheta[x - b]) - 
   P DiracDelta[x - c] + Rbb DiracDelta[x - L];
ode = v''''[x] == load[x]/EI
c = {v[0] == 0, v[L] == 0, v''[0] == 0, v''[L] == 0};

int1 = Integrate[ode[[2]], x, 
         Assumptions -> {x, x1, x2, x3, x4} \[Element] Reals] + x1

int2 = Integrate[int1, x, 
         Assumptions -> {x, x1, x2, x3, x4} \[Element] Reals] + x2

int3 = Integrate[int2, x, 
         Assumptions -> {x, x1, x2, x3, x4} \[Element] Reals] + x3 

int4[x_, x1_, x2_, x3_, x4_, L_, P_, EI_] = 
        Integrate[int3, x, 
       Assumptions -> {x, x1, x2, x3, x4} \[Element] Reals] + x4

With that many HeaviSideTheata, Solve works only, if you insert definite L. You can take any value. Means, at that at a definite number of points function int4 is not defined. But i think, that's no problem.

sol = Flatten@
   Solve[{int4[0, x1, x2, x3, x4, L, P, EI] == 0, 
      int4[L, x1, x2, x3, x4, L, P, EI] == 0, 
      Derivative[2, 0, 0, 0, 0, 0, 0, 0][int4][0, x1, x2, x3, x4, L, 
        P, EI] == 0, 
      Derivative[2, 0, 0, 0, 0, 0, 0, 0][int4][L, x1, x2, x3, x4, L, 
        P, EI] == 0} /. L -> EulerGamma, {x1, x2, x3, x4}, Reals] /. 
  EulerGamma -> L

(*   {x1 -> -(P/EI), x2 -> 0, x3 -> -((335 L^2 P)/(6144 EI)), x4 -> 0}   *)

vsol[L_, P_, EI_][x_] = 
   Evaluate[int4[x, x1, x2, x3, x4, L, P, EI] /. sol]

(*   (1/(6144 EI L))(-335 L^3 P x + 416 L P x^3 HeavisideTheta[x] - 
  864 L P (L - x)^3 HeavisideTheta[-L + x] + 
  16 L^4 P HeavisideTheta[-(L/2) + x] - 
  128 L^3 P x HeavisideTheta[-(L/2) + x] + 
  384 L^2 P x^2 HeavisideTheta[-(L/2) + x] - 
  512 L P x^3 HeavisideTheta[-(L/2) + x] + 
  256 P x^4 HeavisideTheta[-(L/2) + x] - 
  L^4 P HeavisideTheta[-(L/4) + x] + 
  16 L^3 P x HeavisideTheta[-(L/4) + x] - 
  96 L^2 P x^2 HeavisideTheta[-(L/4) + x] + 
  256 L P x^3 HeavisideTheta[-(L/4) + x] - 
  256 P x^4 HeavisideTheta[-(L/4) + x] + 
  432 L^4 P HeavisideTheta[-3 L + 4 x] - 
  1728 L^3 P x HeavisideTheta[-3 L + 4 x] + 
  2304 L^2 P x^2 HeavisideTheta[-3 L + 4 x] - 
  1024 L P x^3 HeavisideTheta[-3 L + 4 x])   *)
Manipulate[
 Plot[{vsol[L, P, EI][x], vsol[L, P, EI]'[x], vsol[L, P, EI]''[x], 
   vsol[L, P, EI]''''[x]}, {x, -5, 5}, PlotRange -> 5, 
  PlotStyle -> {Red, Green, Blue, {Thick, Magenta}}], {{L, 3}, -4, 
  4}, {{P, 2}, -4, 4}, {{EI, 1}, -4, 4}]

For the test, whether lhs - rhs of ode is zero, you again get a result with a few points excluded.

(Subtract @@ ode) /. v -> vsol[L, P, EI] /. L -> EulerGamma // 
   Expand // FullSimplify

(*   0   *)

Plot[Evaluate[((Subtract @@ ode) /. v -> vsol[L, P, EI] // Expand) /. 
   L -> EulerGamma], {x, -5, 5}, PlotStyle -> Thick]