Solve a system of second order PDE

[Woppe]

I need to solve the following system of differential equations:
[tex]\frac{\partial^2 y}{\partial t^2} + A\frac{\partial y}{\partial t} - B \frac{\partial^2 y}{\partial z^2} = Cq[/tex]

[tex]
\frac{\partial^2 q}{\partial t^2} + D\frac{\partial q}{\partial t} + q = E\frac{\partial^2 y}{\partial t^2}
[/tex]

A,B,C,D and E are constants.
Is it possible to solve this using Matlab? In that case, how is it done? I have looked into the function pdepe but I could not figure out how to use it in this case.

 

[kosovtsov]

Present Computer Algebra Systems are weak enough in automatic PDEs solving. At the same time they provide unprecedented possibilities to carry out exhausting calculations.

What about your problem. First of all you need split your system. Finding [tex]q[/tex] from first DE and substituting it into second DE gives the equivalent PDE system, where the second PDE has only one dependent variable. It’ll be linear fourth order in [tex]t[/tex] and second order in [tex]z[/tex].

Seek particular solution of the second PDE in the following form

[tex]y(t,z) = F(\tau,z)\exp(\tau t)[/tex]

Substituting it, you obtain the second order DE for F (in fact it is ODE) which is solvable. Then (the general) solution to second PDE looks like

[tex]y(t,z) = \int_{- \infty}^\infty F(\tau,z)\exp(\tau t)d \tau[/tex]

where, of course, [tex]F[/tex] is the solution of above DE.

Substitution of it in the first DE gives the solution for [tex]q[/tex].

 

[Woppe]

I understand how you create the new DE, but then I don't understand how you denote your expression.

Am I supposed to find F and tau? What substitution is required to obtain the second expression?

Thanks in advance.

 

[elibj123]

kosovtsov actually performed some kind of Fourier/Laplace transform.

What is suggested is that you will transform the time-domain (t variable) to a frequency domain (w variable) so your function y(t,z) is represented:

[tex]y(t,z)=\int^{\infty}_{-\infty}F(w,z)e^{iwt}dw[/tex]

If you're not quite familiar with Fourier transform, I will suggest you reading into it - but generally it means that your function y is composed of infinite harmonics with different frequency (exp(iwt) ) and different amplitudes which depend on the frequency (F(w)).
The function F(w) (in our case its F(w,z)) is the Fourier transform of y(t,z).

What is comfortable about Fourier transform, is that differentiation in the time domain, becomes multiplication in the frequency domain. So:

[tex](\partial_{t})^{n}y(t,z)=(iw)^{n}F(w,z)[/tex]
[tex](\partial_{z})y(t,z)=\partial_{z}F(w,z)[/tex]

So you reduce your PDE to an ODE-like equation (because you can regard w as a parameter, not a variable you differentiate wrt).

So go over your new PDE, and whenever you see y swap it with F, and whenever you see a t-derivative swap it with a powered (iw) and then you will have a simpler equation for F. Of course this doesn't mean you will have an explicit solution to y (it'll be expressed as an integral) but this can be approximated with numerical methods.

 

[Woppe]

I have been talking to a professor in numerical calculations at my university, and I was told this equation is difficult to solve. I was recommended to look into second order wave equations, solved by finite differences. I have still not come closer to solving the problem though.

How can I find F in the equations you have written?

 

[kosovtsov]

The equivalent splited system is as follows

[tex]
q= (\frac {\partial^2 y}{\partial t^2}+A\frac {\partial y}{\partial t}-B\frac {\partial ^2 y}{\partial z^2})/C
[/tex]

[tex] \frac {\partial ^4y}{\partial t^4}= B\frac {\partial ^4y}{\partial t^2 \partial z^2}-A\frac {\partial ^3y}{\partial t^3}-AD\frac {\partial^2y}{\partial t^2}+\frac {\partial ^2y}{\partial z^2}-A\frac {\partial y}{\partial t}-\frac {\partial^3y}{\partial t^3}+CE\frac {\partial^2y}{\partial t^2} -\frac {\partial^2y}{\partial t^2}+B\frac {\partial^2y}{\partial z^2} [/tex]

Seeking a partial solution in form

[tex] y(t,z) = F(\tau,z)\exp(\tau t) [/tex]

we get the following DE for [tex] F(\tau,z)[/tex]

[tex] -B(\tau^2+D\tau+1)\frac {\partial ^2F}{\partial z^2}+\tau(\tau^3-EC\tau+D\tau^2+\tau+A\tau^2+AD\tau+A)F[/tex]

In fact it is the oscillator equation (the scond order ODE with constant coefficients, as coefficients do not depend on [tex]z[/tex]. Its solution is

[tex] F(\tau,z) = F1(\tau)\exp(\sqrt{\frac{\tau(\tau^3-EC\tau+D\tau^2+\tau+A\tau^2+AD\tau+A)}{B(\tau^2+D\tau+1)}}z)+\\ F2(\tau)\exp(-\sqrt{\frac{\tau(\tau^3-EC\tau+D\tau^2+\tau+A\tau^2+AD\tau+A)}{B(\tau^2+D\tau+1)}}z) [/tex]

So the general solution for [tex]y[/tex]

[tex] y(t,z) = \int_{-\infty}^\infty F(\tau,z)\exp(\tau t) d\tau\,,[/tex]

where [tex] F(\tau,z)[/tex] is in above form and [tex] F1(\tau), F2(\tau)[/tex] are arbitrary functions.

 

[Woppe]

Thanks for that!

To calculate the integral, I wrote a program in Matlab.
I assign the constans A, B, C, D, E values between 0.1 and 1.0, and use the following command:

y = quadgk(@(tau)FsfaTauZ(tau,z).*exp(tau.*t),-Inf,Inf)

Where FsfaTauZ contains the the function F as you wrote it, with F1 and F2 being constants equal to 1.
Just for the sake of simplicity, also t and z are set to 1.
Problem is, the value returned by quadgk approaches infitity and has a small imaginary component as well. That's not reasonable, since the value of y describes the deflection of a line exposed to a wind load.