Help with difficult differential equation

[hushish]

Hi,

I need help with a difficult differential equation. The following is the differential equation for a beam on a varying-stiffness elastic foundation. Essentially the beam sits on a cantilever, whose stiffness varies to the third power with length. Does anybody know how to solve such an equation as a closed form solution (if possible)?

E*I*y(x)''''+k*y(x)/x^3=0

where k is a constant.

 

[tiny-tim]

welcome to pf!

hi hushish! welcome to pf!

(try using the X² button just above the Reply box )

have you tried y = a power of x ?

 

[chiro]

Hi,

I need help with a difficult differential equation. The following is the differential equation for a beam on a varying-stiffness elastic foundation. Essentially the beam sits on a cantilever, whose stiffness varies to the third power with length. Does anybody know how to solve such an equation as a closed form solution (if possible)?

E*I*y(x)''''+k*y(x)/x^3=0

where k is a constant.

Hey hushish and welcome to the forums.

Are you familiar with integral transforms applied to linear differential equations? Have you ever come across the Laplace transform or the Fourier Transform?

 

[hushish]

I haven't done Laplace transforms since University. I'm open to anything that will yield a solution, but I'll need to be guided through it...

 

[LCKurtz]

If you lump the constants together and multiply the equation by $x^4$ it is of the form$$
x^4y'''' +cxy=0$$That is an Euler - Cauchy equation. See

http://en.wikipedia.org/wiki/Cauchy–Euler_equation

which discusses various ways to solve it. One of the most direct is to follow Tiny Tim's suggestion of looking for a solution of the form $y = x^p$. You don't want to mess with LaPlace transforms for this one.

 

[tiny-tim]

You don't want to mess with LaPlace …

'nuff said!

 

[hushish]

Thanks for the advice and the link. I am now left with the following equation:

m⁴-6m³+11m²-6m+kx=0

Any suggestions for the roots of the equation? Some of them will be complex...a little too much maths for my engineering brain.

 

[LCKurtz]

Thanks for the advice and the link. I am now left with the following equation:

m⁴-6m³+11m²-6m+kx=0

Any suggestions for the roots of the equation? Some of them will be complex...a little too much maths for my engineering brain.

You mean $m^4-6m^3+11m^2-6m+k$ with no $x$. I would suggest a mathematics program such as Maple of Mathematica. Maple seems to be able to solve it explicitly. If you don't have access to such a program, tell me what $k$ is and I will give you the roots. Or are you hoping for a general formula containing $k$ as a parameter?

 

[pasmith]

If you lump the constants together and multiply the equation by $x^4$ it is of the form$$
x^4y'''' +cxy=0$$That is an Euler - Cauchy equation.

Except it isn't, because $xy$ is not of the form $x^ny^{(n)}$.
This means that substituting $y = x^m$ doesn't work because it leads to:

m⁴-6m³+11m²-6m+kx=0

Any suggestions for the roots of the equation? Some of them will be complex...a little too much maths for my engineering brain.

Don't bother: If $y = x^m$ were actually a solution, you would have ended up with a fourth-order polynomial for $m$ whose coefficients are independent of $x$. (As an aside, a pair of complex conjugate roots $m = p \pm iq$ would result in real-valued solutions of the form $x^p \cos(q\ln(x))$ and $x^p \sin(q\ln(x))$).

I suspect the only method of solving this one (aside from doing it numerically, which would be my first choice here) is to use a variant of Frobenius' Method: pose a power series solution of the form
$$y(x) = x^k \sum_{n=0}^{\infty} a_nx^n$$
and choose $k$ so that $a_0$ is not required to be zero. Substituting this into
the ODE then gives the condition
$$a_0 k(k-1)(k-2)(k-3) = 0$$
Unfortunately the roots of $k(k-1)(k-2)(k-3) = 0$ differ by integers, so we only get one linearly independent solution from this; and for $a_1$ to be finite we must take $k = 3$. We then have the recurrence relation
$$a_{n+1} = {{-c a_n} \over {(n+4)(n+3)(n+2)(n+1)}}<br /> = -c a_n {{n!} \over {(n+4)!}}$$
which gives
$$a_n = a_0 (-1)^n c^n \prod_{r=0}^{n-1} {{r!} \over {(r+4)!}}$$
The resulting power series has an infinite radius of convergence: for $n \geq 4$, we have
$$|a_0| |c|^n |x|^{n+3} \prod_{r=0}^{n-1} {{r!} \over {(r+4)!}} = <br /> |a_0| |x|^3 |cx|^n {{2!3!} \over {n!(n+1)!(n+2)!(n+3)!}} <<br /> |a_0| |x|^3 2!3! {{|cx|^n} \over {n!}}$$
Thus, for $n \geq 4$, the absolute value of the $n$th term of the series for $y(x)$ is strictly less than a constant ($12|a_0| |x|^3$) times the $n$th term of the series for $e^{|cx|}$, which converges for all finite $|cx|$. Since the behaviour of the first three terms does not affect convergence, the series for $y(x)$ also converges for all finite $|cx|$.

That gives one solution; there are three others to find.

 

[LCKurtz]

Except it isn't, because $xy$ is not of the form $x^ny^{(n)}$.

You're absolutely correct. Big woops on my part.

 

[hushish]

Thanks a mil guys. That seems like the way to go. I have a few understanding questions on the maths, and if I come across any problems I'll post it here.

 

[pasmith]

On further reflection, I think the only ways of getting all four linearly independent solutions are (1) an integral transform method (possibly Laplace) which hopefully will yield four separate contours on the complex plane each of which will contribute a solution, or (2) solve it by fiat: give names to the following four special cases and declare that the solution is a linear combination of them: (a) $y(1) = 1$, $y'(1) = y''(1) = y'''(1) = 0$, (b) $y(1) = y''(1) = y'''(1) = 0$, $y'(1) = 1$ , (c) $y(1) = y'(1) = y'''(1) = 0$, $y''(1) = 1$, and (d) $y(1) = y'(1) = y''(1) = 0$, $y'''(1) = 1$. Obviously, if you're interested in the actual behaviour or values of these functions you'll either have to use the first method or calculate them by numerical solution of the ODE. But frankly, if you're reduced to solving it numerically you may as well just solve it for the particular initial or boundary conditions you have. (Be aware that $x = 0$ is a singular point of the ODE, so one or more of these might blow up there. That's why I'm applying conditions at $x = 1$, where all four should be well-behaved.)

 

[hushish]

Thanks pasmith. Unfortunately the maths that you describe is beyond me. I'm just a simple engineer, and so I have never manipulated the maths in the way you describe. Is there any way that we can work together on the solution somehow-perhaps you can guide me in the steps. In return, I can source you in any documentation that comes out as a result of the research...

 

[chiro]

Thanks pasmith. Unfortunately the maths that you describe is beyond me. I'm just a simple engineer, and so I have never manipulated the maths in the way you describe. Is there any way that we can work together on the solution somehow-perhaps you can guide me in the steps. In return, I can source you in any documentation that comes out as a result of the research...

For your equation: E*I*y(x)''''+k*y(x)/x^3=0, you can transform this to x^3*E*I*y(x)''''+k*y(x)=0, where you will have to not use x=0 (unless your initial condition y(0) = 0 holds) and like hushish said, use Laplace techniques.

You say you are an engineer, so I'm surprised you have not come across the Laplace transform technique before.

If you don't need to prove it rigorously then use something like Wolfram Alpha, but if you want to prove it an understand it, then you will really need to learn the Laplace Transform and the identities for the forward transformation and for the backward (inverse) transformation.

You can find everything you need to do for forward transformation by looking at a decent table of Laplace transformations and also by the definition of the Laplace transformation by it's integral representation and then deriving results using normal calculus techniques that usually involve integration by parts.

 

[Dickfore]

If you lump the constants together and multiply the equation by $x^4$ it is of the form$$
x^4y'''' +cxy=0$$That is an Euler - Cauchy equation. See

http://en.wikipedia.org/wiki/Cauchy–Euler_equation

which discusses various ways to solve it. One of the most direct is to follow Tiny Tim's suggestion of looking for a solution of the form $y = x^p$. You don't want to mess with LaPlace transforms for this one.

No, it is not. An Euler equation would be:
$$x^4 y^{\mathrm{IV}} + c \, x \, y' = 0$$
or
$$x^4 \, y^{\mathrm{IV}} + c \, y = 0$$

 

[LCKurtz]

No, it is not. An Euler equation would be:
$$x^4 y^{\mathrm{IV}} + c \, x \, y' = 0$$
or
$$x^4 \, y^{\mathrm{IV}} + c \, y = 0$$

As was noted earlier and I have already acknowledged a couple of days ago in post #10. Is there some reason to point it out again?

 

[pasmith]

For your equation: E*I*y(x)''''+k*y(x)/x^3=0, you can transform this to x^3*E*I*y(x)''''+k*y(x)=0, where you will have to not use x=0 (unless your initial condition y(0) = 0 holds) and like hushish said, use Laplace techniques.

You say you are an engineer, so I'm surprised you have not come across the Laplace transform technique before.

If you don't need to prove it rigorously then use something like Wolfram Alpha, but if you want to prove it an understand it, then you will really need to learn the Laplace Transform and the identities for the forward transformation and for the backward (inverse) transformation.

You can find everything you need to do for forward transformation by looking at a decent table of Laplace transformations and also by the definition of the Laplace transformation by it's integral representation and then deriving results using normal calculus techniques that usually involve integration by parts.

That's not really helpful here: the equation has non-constant coefficients. Indeed Wolfram Alpha can't handle the Laplace transform of $x^3y^{(4)}$: (annoyingly I can't include links, so go to Wolfram alpha, type in "laplace transform of t^3f''''(t)" and see what you get). Nor can Wolfram alpha solve the ODE directly.

What you have to do is to let
$$y(x) = \int_{\gamma} e^{-px} f(p) dp$$
(I don't know why I used e^{-px} as my kernel; e^{px} would have been the obvious choice)
where f(p) is to be determined and $\gamma$ is a contour in the complex p-plane which is to be chosen to our advantage. Substituting this, we end up with
$$[ e^{-px} g(x,p)]_{\gamma} + \int_{\gamma} e^{-px} D(f) dp = 0$$
where g(x,p) is obtained by integration by parts and D is a differential operator with respect to p. The idea is then to solve $D(f) = 0$ for $f(p)$, and then to choose a contour $\gamma$ so that $\gamma$ encloses one or more singularities of $f$ and the boundary term $[ e^{-px} g(x,p)]_{\gamma}$ vanishes (hopefully choosing different contours will give us the four linearly independent solutions of the original ODE).

For the equation in question, $D(f)$ reduces to $p^4u'''(p) + cu(p)$ where $u(p) = p^4 f(p)$ which is a slight improvement (we've reduced it to third order) and still linear but still not straightforward to solve.

Really, I think it would be better to spend time and effort on writing a program which will take initial/boundary conditions and constants and solve the ODE numerically, rather than trying to solve it analytically.

 

[pasmith]

I believe I have found all four solutions.

We have:
$$y^{(4)} + {{ky} \over {x^3}} = 0$$

My first step was to substitute $y = x^3f(x)$ which gave

$$x^3 f^{(4)} + 12x^2 f^{(3)} + 36xf'' + 24f' + kf = 0$$

I then looked for a power series solution for f of the form
$$f(x) = x^{\lambda}\sum_{n=0} a_n(\lambda) x^n$$
where $\lambda$ will be chosen such that $a_0 \neq 0$. Substituting the above into the ODE for f gave the following condition on $\lambda$:

$$\lambda(\lambda^3 + 6\lambda^2 + 47\lambda - 30) = 0$$

The solution $\lambda = 0$ was not a surprise: it corresponds to the series I found before. But the roots of the cubic factor are new, and according to Wolfram alpha are not integers: there's a real solution $\lambda_0 ≈ 0.6$ and a complex conjugate pair $p \pm iq$ with a negative non-integer real part.

The recurrence relation for $a_n(\lambda)$, $n \geq 1$, is
$$a_n(\lambda) = {{-ka_{n-1}(\lambda)} \over {(n+\lambda)((n+\lambda)^3 + 6(n+\lambda)^2 + 47(n+\lambda) - 30)}}$$

For the complex roots, it's not obvious that $a_n(\lambda)$ is real if $a_0(\lambda)$ is real. However we can set $a_n(\lambda) = b_n(\lambda) + ic_n(\lambda)$ for real $b_n(\lambda)$ and $c_n(\lambda)$. And in fact if we take $a_0(p+iq) = A + iB$ and $a_0(p-iq)=A-iB$ for arbitrary real A and B, then the coefficients obtained from $\lambda = p-iq$ will be the complex conjugates of those obtained from $\lambda = p+iq$. So we get the same $b_n$ for each, with $c_n$ differing only by a sign, which we fix by considering $\lambda = p+iq$. We then have $b_0 = A$, $c_0 = B$.

Taking $\lambda = p+iq$ gives
$$f(x) = x^{p+iq} \sum_{n=0}^{\infty} (b_n + ic_n)x^n<br /> = x^p e^{iq\ln x} \sum_{n=0}^{\infty} (b_n + ic_n)x^n<br /> = x^p (\cos(q\ln x) + i\sin(q\ln x)) \sum_{n=0}^{\infty} (b_n + ic_n)x^n$$
whose real and imaginary parts (provided $x > 0$) are
$$\Re(f(x)) = x^p\cos(q\ln x) \sum_{n=0}^{\infty} b_n x^n - <br /> x^p\sin(q\ln x) \sum_{n=0}^{\infty} c_n x^n$$
and
$$\Im(f(x)) = x^p\sin(q\ln x) \sum_{n=0}^{\infty} b_n x^n + <br /> x^p\cos(q\ln x) \sum_{n=0}^{\infty} c_n x^n$$
Taking $\lambda = p-iq$ simply gives the complex conjugate of the above.

So the four real-valued solutions are
$$\lambda = 0: y(x) =x^3f(x) = x^3 \sum_{n=0}^{\infty} a_n(0)x^n$$
$$\lambda = \lambda_0: y(x) =x^3f(x) = x^{3+\lambda_0} \sum_{n=0}^{\infty} a_n(\lambda_0)x^n$$
$$\lambda = p \pm iq: y(x) =x^3f(x) = x^{p+3}\cos(q\ln x) \sum_{n=0}^{\infty} b_n x^n - <br /> x^{p+3}\sin(q\ln x) \sum_{n=0}^{\infty} c_n x^n$$
$$\lambda = p \pm iq: y(x) =x^3f(x) = x^{p+3}\sin(q\ln x) \sum_{n=0}^{\infty} b_n x^n + x^{p+3}\cos(q\ln x) \sum_{n=0}^{\infty} c_n x^n$$
where the last two are valid only for $x > 0$.

I don't propose to solve the recurrence relations, but I can prove convergence of the power series by the ratio test:

$$\lim_{n \to \infty} \left| {{a_n(\lambda) x^n} \over {a_{n-1}(\lambda)x^{n-1}}} \right|<br /> = \lim_{n \to \infty} \left| {{kx} \over {(n+\lambda)((n+\lambda)^3 + 6(n+\lambda)^2 + 47(n+\lambda) - 30)}}\right| = 0$$

It's satisfying to have found all four solutions analytically, but I still think that, in practice, numerical solution is the way to go.

 

[Dickfore]

The point $x = 0$ is singular for the equation:
$$x^3 y^{IV} = -c y$$
because it implies
$$y(0) = y'(0) = y''(0) = 0$$
as can easily be verified by differentiating twice and substituting $x = 0$. This means we cannot give arbitrary initial conditions to find the 4 conditions in a general solution. We can only choose one arbitrary constant corresponding to the freedom in choosing $y'''(0)$. Thus, we would find a singular solution.

Taking nth derivative and setting $x = 0$, we get:
$$6 \, \left( \begin{array}{c}n \\ 3\end{array} \right) \, y^{(n + 1)}(0) = -c y^{(n)}(0), \ n \ge 3$$

This first order recursion allows us to write down a Taylor series solution. Take $c_{n} = y^{(n + 3)}/(n + 3)!$. We get:
$$6 \frac{(n + 3)!}{3! \, n!} (n + 4)! c_{n + 1} = -c (n + 3)! c_{n}, \ n \ ge 0$$

$$c_{n+1} = -c \frac{c_{n}}{(n + 1) (n + 2) (n + 3)}, \ n \ge 0$$
The solution of this recursion is:
$$c_{n} = c_{0} \frac{0! \, 1! \, 2! \, (-c)^{n}}{n! \, (n + 1)! \, (n + 2)!}$$
The series solution is
$$y(x) = K \, x^3 \, f(c \, x)$$
where:
$$f(x) = \sum_{n = 0}^{\infty}{\frac{(-x)^{n}}{n! \, (n + 1)! \, (n + 2)!}}$$

 

[Dickfore]

Here are 3 plots of $x^3 \, f(x)$. The sum should converge rapidly. I calculated the sum of the first 201 terms by Mathematica. The function is oscillatory with ever increasing amplitude. Notice the change of vertical scale:

 

[pasmith]

The point $x = 0$ is singular for the equation:
$$x^3 y^{IV} = -c y$$
because it implies
$$y(0) = y'(0) = y''(0) = 0$$
as can easily be verified by differentiating twice and substituting $x = 0$. This means we cannot give arbitrary initial conditions to find the 4 conditions in a general solution.

This is a fourth-order linear ODE. There are four linearly-independent solutions. Zero is a singular point, so we should expect some of these to have singularities there. And indeed two of them do: the ones which behave like $x^p\cos(q\ln x)$ and $x^p\sin(q\ln x)$ oscillate with increasing frequency and amplitude as $x \to 0$. But the other two are well-behaved and vanish. Thus at the origin we can only impose one condition on $y'$ (plus the second condition "bounded at origin"). On the other hand, if the origin isn't in the domain in which we want the solution to make physical sense, all four solutions are available. (Compare, by way of example, Bessel's equation
$$x^2 y'' + xy' + x^2y = 0$$
which has a singular point at $x = 0$ and two linearly dependent solutions: $J_0(x)$, which is finite at the origin ($J_0(0) = 1$), and $Y_0(x)$, which has a singularity.)

The series solution is
$$y(x) = K \, x^3 \, f(c \, x)$$
where:
$$f(x) = \sum_{n = 0}^{\infty}{\frac{(-x)^{n}}{n! \, (n + 1)! \, (n + 2)!}}$$

I don't think so; using that I get (with K = 1):
$$y^{(4)} = \sum_{n=0}^{\infty} {{(-c)^{n+1}(n+4)x^{n}} \over {n!(n+1)!(n+2)!}}$$
whereas
$$-{cy \over x^3} = \sum_{n = 0}^{\infty}{\frac{(-c)^{n+1}x^n}{n! (n + 1)! (n + 2)!}}$$
and the two are not equal. I derived the obvious series solution in the post in which I originally pointed out that we were not dealing with an Euler equation:
$$y(x) = x^3\sum_{n=0}^{\infty} a_nx^n = a_0x^3\sum_{n=0}^{\infty} \frac{2!3!(-1)^nc^nx^n}{n!(n+1)!(n+2)!(n+3)!}$$

 

[pasmith]

The function is oscillatory with ever increasing amplitude.

You can't draw that conclusion from a truncated series: what you have is a polynomial of order 203 with a triple repeated root at the origin. It's going to appear oscillatory (there are 200 roots aside from the origin; I'd expect a good number of them to be real and positive) and tend to infinity as x tends to infinity. The same can be said of the truncation of the series for exp(-x) after 201 terms: it's an order 200 polynomial, so one expects it to appear oscillatory and to tend to infinity as x tends to infinity; but we know that exp(-x) tends to zero as x tends to infinity and has no real roots.

If you want to look at large $x$, you need to substitute $x = 1/z$ in the ODE and look at the limit $z \to 0$.

 

[Dickfore]

You can't draw that conclusion from a truncated series...

Yes, you can. For example:

As you increase the number of terms in the Taylor series for $\sin x$, the domain where the truncated series approximates the true function increases, eventually encompassing the zeros of the function.

The oscillatory nature of the function is seen from the alternating sign of successive terms in the series (see Bessel function, for another example).

 

[Dickfore]

You are right, the recursion relation for the coefficients:
$$c_{n} \equiv \frac{y^{(n + 3)}}{(n + 3)!}$$
should read:
$$\frac{(n + 3)!}{n!} \, (n + 4)! \, c_{n + 1} = -c \, (n + 3)! \, c_{n}$$
$$c_{n + 1} = \frac{(-c) \, c_{n}}{(n + 1)(n + 2) (n +3) (n + 4)}$$
which has a solution:
$$c_{n} = c_{0} \, \frac{0! \, 1! \, 2! \, 3! \, (-c)^{n}}{n! \, (n + 1)! \, (n + 2)! \, (n + 3)!}, \ n \ge 0$$
This gives the following analytic solution of the given equation:
$$y = K \, x^3 \, f(c \, x)$$
where:
$$f(x) = \sum_{n = 0}^{\infty}{ \frac{(-x)^{n}}{n! \, (n + 1)! \, (n + 2)! \, (n + 3)!}}$$
This series converges even more rapidly.

Here is the partial sum of the first 201, and the first 501 terms, respectively for f(x):

​

As you can see, they are indistinguishable in the much larger region $x \in [0, 200]$. Also, you can see the function still changes sign.

 

[hushish]

Thanks guys. Still above my head, but I'll try solve it using the info you gave me.

 

[haruspex]

hushish,
Pardon me for questioning this, but are you quite sure your initial ODE is correct? It bothers me that it's not automatically dimensionally consistent. That's relatively unusual in mechanics, no?

 

[hushish]

Hi haruspex. It is dimensionally consistent in that within the k constant is the product of the Young's Modulus E (N/mm²) and the inertia I (mm⁴). If you divide this by the cube of the distance we are left with N/mm, and as the foundation constant is per mm width, the dimensions of k is N/mm²-dimensionally consistent.