# Differential equation with a Heavisidetheta function running for ever

So I'm trying to get the solution of the following differential equation:

DSolve[{
y''[x] - a^2 y''''[x] ==
-(c + (T - c) x/(L/2) - (2 (T - c) (x - L/2)/(L/2)) HeavisideTheta[x - L/2])/(GI),
y[0] == 0, y''[0] == 0, y[L] == 0, y''[L] == 0}, y, x]


Mathematica keeps running without getting any solution.

Anyone knows what may be the problem? Is the syntax right?

Mathematica might need some help with this one.

One of the problems is surely that you are asking for y - try first to get a solution for u(x)=y''(x), this makes things a bit easier. You can always integrate twice at the end.

In this case, it might be an idea to split up the 'source term' (i.e. the right hand side). As the differential equation is linear, you can always add the solutions.

So I get

deq1 = {u[x] - a^2 u''[x] == -(c + (T - c) x/(L/2))/GI}
deq2 = {u[x] - a^2 u''[x] == -(-2*(T - c) (x - L/2)/(L/2)*HeavisideTheta[x - L/2])/GI}


and then

sol1 = First@DSolve[deq1, u[x], x]
sol2 = First@DSolve[deq2, u[x], x]


Put the solutions together again:

uu = ((u[x] /. sol1) + (u[x] /. sol2))


Now solve for your boundary conditions y''[0]=y''[L]=0 (or u[0]=u[L]=0):

Solve[ { (uu /. x -> 0 ) == 0, (uu /. x -> L) == 0}, {C[1], C[2]}]
uu /. First@%
uu = Simplify[%]


The result is longish, but now you have your solution for u(x) (or y''(x))

The next steps are to integrate twice over x (don't forget the integration constants)

y = Integrate[uu,x,x] + C[1] x + C[2]


and use the boundary conditions for y(0) and y(L) to get rid of the integration constants:

y = Simplify[y /. First@Solve[{(y /. x -> 0) == 0, (y /. x -> L) == 0} , {C[1], C[2]}]]


Still not a very nice result, but there you are.

If Mathematica can't decide if argument to HeavisideTheta is positive or negative, this remains unevaluated. Type HeavisideTheta[x - L/2] on your computer and see the result. So, the result will contain this even if it can solve it. I think it just can't find a particular solution because of this.

One way to help M, is to break it to 2 ODE to avoid the explicit use of unit step in the RHS.

One ODE is solved for first half, and the second ODE for the second half. You use the conditions at the end of the first half from the first solution as initial conditions for the second half.

Like this:

solFirstHalf = y[x] /. First@DSolve[{
y''[x] - a^2 y''''[x] == -(c + (T - c) x/(L/2))/(GI),
y[0] == 0, y''[0] == 0, y[L] == 0, y''[L] == 0},
y[x], x];


Now solve on the second half. Now we add the stuff that had the shifted unit step, but we throw away the unit step function as not needed

yAtHalf = solFirstHalf /. x -> L/2;
yDerAtHalf = D[solFirstHalf, {x, 2}] /. x -> L/2
solSecondtHalf = y[x] /. First@DSolve[{
y''[x] - a^2 y''''[x] == -(c + (T - c) x/(L/2) - (2*(T - c) (x - L/2)/(L/2)))/(GI),
y[L/2] == yAtHalf, y''[L/2] == yDerAtHalf, y[L] == 0, y''[L] == 0},
y[x], x]


myFinalSol[x_, L_] := Piecewise[{{solFirstHalf, x <= L/2}, {solSecondtHalf, True}}];
vals = {L -> 1, c -> 2, T -> 3, GI -> 4, a -> 5};
Plot[myFinalSol[x, 1] /. vals, {x, 0, 1}]


To verify, compare to numerical solution (now add the unitstep, it is ok, this is numerical)

L = 1;    c = 2;    T = 3;    GI = 4;    a = 5;
sol = y /. First@NDSolve[{y''[x] -
a^2 y''''[x] == -(c + (T - c) x/(L/2) -
(2*(T - c) (x - L/2)/(L/2))*HeavisideTheta[x - L/2])/(GI),
y[0] == 0, y''[0] == 0, y[L] == 0, y''[L] == 0}, y, x];

Plot[sol[x], {x, 0, 1}]


Here is your analytical solution btw. The first valid for first half of the beam x<L/2 and the second is for second half

• You've made a small mistake. y[L] == 0 and y''[L] == 0 can't be used in the first half as boundarys. If you Show the numeric and symbolic result, you'll see they're different :) – xzczd Jan 22 '15 at 8:55

A workaround based on LaplaceTransform:

LaplaceTransform[
y''[x] - a^2 y''''[x] ==
-(c + (T - c) x/(L/2) - (2 (T - c) (x - L/2)/(L/2)) HeavisideTheta[x - L/2])/GI, x, s] /.
Rule @@@ {y[0] == 0, y''[0] == 0};
Solve[%, LaplaceTransform[y[x], x, s]][[1, 1, 2]];
yhalf[x_] = InverseLaplaceTransform[%, s, x];
sol[x_] = yhalf[x] /. Solve[{y[L] == 0, y''[L] == 0} /. y -> yhalf, {y'[0], y'''[0]}][[1]];
vals = {L -> 1, c -> 2, T -> 3, GI -> 4, a -> 5};
(* Evaluate is added in Plot for speed *)
Plot[sol[x] /. vals // Evaluate, {x, 0, 1}]