0
$\begingroup$

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].

$\endgroup$
1
  • $\begingroup$ You have numerous basic coding errors here. As far as I see, every one of your assignments is incorrect. (Should be f[x_] := ..., for example.). You should also avoid variables starting with upper-case letters as they conflict with Mathematica's function names. And more. $\endgroup$ Mar 31 at 7:30

1 Answer 1

0
$\begingroup$

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}]

enter image description here

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]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.