Solving a PDE : 2 order in time, 4 order in space, mixed derivatives

[vatant]

I have a problem that I tried to solve using MAPLE but I guess wasnt doing the right thing.

\frac {\partial ^{2} \delta}{\partial t^{2}}+ S*(\frac {\partial^{2}\delta}{\partial \eta^{2}}+M*\frac {\partial^{4}\delta}{\partial \eta^{4}})-G* \frac {\partial ^{3} \delta}{\partial \eta ^{2} \partial {t}}=0

A, S, M, G are constants


I have BC
\delta (0,t)=0


\frac{\partial \delta}{\partial \eta} (0,t)=A


Further,

\frac {\partial ^{2} \delta}{\partial z^{2}} (0,t)=0

\frac {\partial ^{3} \delta}{\partial z^{3}} (0,t)=0

\frac {\partial ^{4} \delta}{\partial z^{4}} (0,t)=0

How can i solve this problem> Can maple or MATLAB help?

Any other BC can be provided.

Can somebody help?


Thanks


Ivde

 

[LeBrad]

Can maple or MATLAB help?

Sure, assuming eta is space, discretize it spatially and integrate a big system of equations through time. What is your domain?

 

[saltydog]

I used Mathematica to solve the following IBVP:

\text{DE:}\quad \frac {\partial ^{2} u}{\partial t^{2}}-\frac {\partial^{2}u}{\partial x^{2}}+\frac {\partial^{4}u}{\partial x^4}}-\frac {\partial ^{3} u}{\partial x^{2} \partial {t}}=0\quad 0\le x\le 1

\text{BC:}\quad u(t,0)=u(t,1)=0

\text{BC:}\quad u_x(t,0)=\pi

\text{BC:}\quad u_x(t,0)=-\pi

\text{IC:}\quad u(0,x)=Sin(\pi x)

\text{IC:}\quad u_t(0,x)=0

This is the Mathematica code:

sol = NDSolve[{D[u[t, x], {
    t, 2}] - D[u[t, x], {x, 2}] + D[u[t, x], {x, 
          4}] - D[u[t, x], {t, 1}, {x, 2}] == 0,
      
      (* boundary conditions *)
      u[t, 0] == 0,
      u[t, 1] == 0,
      Derivative[0, 1][u][t, 0] == PI,
      Derivative[0, 1][u][t, 1] == -PI,
      
      (* initial conditions *)
      u[0, x] == Sin[PI x],
      Derivative[1, 0][u][0, x] == 0
      
      },
    u, {t, 0, 1}, {x, 0, 1}, PrecisionGoal -> 2]
Plot3D[Evaluate[u[t, x] /. sol[[1, 1]]], {t, 0, 1}, {x, 0, 1}]

A plot of the solution is attached:

 

[saltydog]

So I'd like to know if the solution I got above indeed satisfies the differential equation. So I calculated all the partials from the numeric solution, and back-substituted into the equation. Ideally I should get zero. This is the code I used in mathematica:

[tex]dif[t,x]=\partial_{t,t}u[t,x]-\partial_{x,x}u[t,x]+\partial_{x,x,x,x}u[t,x]-\partial_{t,x,x}u[t,x][/tex]

I then took cross-sections of dif[t,x] and plotted them to see how close it's to zero. For example, I plotted the cross-section at x=0.4 like this:

Plot[dif[t,0.4],{t,0,1}]

The attached plot shows this analysis. It's kinda' close to zero but not good enough for a moon mission.

 

[saltydog]

For the record, I'd like to submit an improvement to the numerical method I used to evaluate the sample PDE I gave above. However I'm not exactly familiar with the method Mathematica used to achieve this improvement. Here's the code:

sol = NDSolve[{D[u[t, x], {
    t, 2}] - D[u[t, x], {x, 2}] + D[u[t, x], {x, 
          4}] -  D[u[t, x], {t, 1}, {x, 2}] == 0,
      
      (* boundary conditions *)
      u[t, 0] == u[t, 1] == 0,
      Derivative[0, 1][u][t, 0] == PI, Derivative[0, 1][u][t, 1] == -PI ,
      
      (* initial conditions *)
      u[0, x] == Sin[  PI  x],
      Derivative[1, 0][u][0, x] == 0
      },
    u, {t, 0, 1}, {x, 0, 1}, PrecisionGoal -> 10, MaxSteps ->
         Infinity, MaxStepSize -> 0.01, Method -> ImplicitRungeKutta]
u[t_, x_] = Evaluate[u[t, x] /. sol[[1, 1]]]
Plot3D[u[t, x], {t, 0, 1}, {x, 0, 1}, PlotPoints -> 50]

Note the method "ImplicitRungeKutta" which Mathematica defines as "families of arbitrary-order implicit Runge-Kutta methods".

I also improved the results by specifying a step size of 0.01. A plot of dif[t,x] is shown below as well as a cross-section at x=0.4. It's now accurate to 3 decimal places, still not good enough for a moon mission but apparently Mathematica can only go to 10,000 data points. If we had to, we (couldn't we?) could write the algorithm from scratch and use 10 million points to achieve what I would consider a minimum for the mission: 10^-7.

 

[vatant]

Mathematica problem

Hi:
I used Mathematica 5.0 to evaluate the sample code that you had written and came across whole bunch of errors..

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.
...
General::stop: Further output of ReplaceAll::reps will be suppressed during \
this calculation
...

etc..


I am not sure if this requires special cell evaluation.

 

[saltydog]

I'm sorry. I input Pi as PI to get the cut-and-paste to work from Mathematica. Here's the code:

sol = NDSolve[{D[u[t, x], {
    t, 2}] - D[u[t, x], {x, 2}] + D[u[t, x], {x, 
          4}] -  D[u[t, x], {t, 1}, {x, 2}] == 0,
      
      (* boundary conditions *)
      u[t, 0] == u[t, 1] == 0,
      Derivative[0, 1][u][t, 0] == Pi, Derivative[0, 1][u][t, 1] == -Pi ,
      
      (* initial conditions *)
      u[0, x] == Sin[  Pi  x],
      Derivative[1, 0][u][0, x] == 0
          PI
      },
    u, {t, 0, 1}, {x, 0, 1}, PrecisionGoal -> 
        10, MaxSteps -> Infinity, MaxStepSize -> 0.01, Method -> \
ImplicitRungeKutta]
u[t_, x_] = Evaluate[u[t, x] /. sol[[1, 1]]]
Plot3D[u[t, x], {t, 0, 1}, {x, 0, 1}, PlotPoints -> 50]

You know of course you can use a real pi symbol in Mathematica but then that messes up the cut-and-paste. Let me know if you still have problems as I do not wish to cause problems like that.

 

[vatant]

changing Coefficients

Hi:

as you spotted that, pi was the problem.

I tried changing the sign of the second term
\frac {\partial ^{2}u}{\partial x^{2}}
from - to +

I guess the system doesn't acknowlege that ! did u try that sign change? says not a Diff eqn..what do you suggest?

 

[saltydog]

Hi:

as you spotted that, pi was the problem.

I tried changing the sign of the second term
\frac {\partial ^{2}u}{\partial x^{2}}
from - to +

I guess the system doesn't acknowlege that ! did u try that sign change? says not a Diff eqn..what do you suggest?

I have tried various changes to the PDE and NDSolve evaluates it accordingly. I suspect you have some other syntactical error with the code. I bet you've already figured that out.

Also, I never did quite understand your original IBVP up there (the eta and z variables). Have you modified the code to fit your particular problem?

 

[vatant]

Pde

Hi there:

I changed the sign of the second terma and NDsolve is not able to comprehend the problem.

your code with a - sign for the second term worked fine. Just by changing the sign of the second term, things are falling apart..is it becos of chanign the PDE from Ellip to Hyper or the other way!

NDSolve[{D[u[t, x], {
t, 2}] + D[u[t, x], {x, 2}] + D[u[t, x], {x,
4}] - D[u[t, x], {t, 1}, {x, 2}] == 0,

I guess, changing the sign has more complications. I watched the problem closely and infact your eqn differs from mine with taht sign and for the preliminary work ..i tried to change the sign of ur code and things blew up.

wud be great to know your thoughts on this.

 

[saltydog]

I have 5.2 and it seems to work. Here's the code with the positive value:

NDSolve[{D[u[t, x], {t, 2}] + D[u[t, x], {x, 2}] +
 D[u[t, x], {x, 4}] -  D[u[t, x], {t, 1}, {x, 2}] == 0,
    
    (* boundary conditions *)
    u[t, 0] == u[t, 1] == 0,
    Derivative[0, 1][u][t, 0] == Pi, Derivative[0, 1][u][t, 1] == -Pi ,
    
    (* initial conditions *)
    u[0, x] == Sin[  Pi  x],
    Derivative[1, 0][u][0, x] == 0
    },
  u, {t, 0, 1}, {x, 0, 1}, PrecisionGoal -> 
        10, MaxSteps -> Infinity, MaxStepSize -> 0.01, 
            Method -> ImplicitRungeKutta]
u[t_, x_] = Evaluate[u[t, x] /. sol[[1, 1]]]
Plot3D[u[t, x], {t, 0, 1}, {x, 0, 1}, PlotPoints -> 50]

Attached is the solution. Notice up there in my previous post I have a stray 'PI' in the code. You took that out right?

 

[vatant]

Mathematica 5.0

I am using version 5.0 of mathematica and looks like it doesn't going to work..you think that's the version problem !?? I should try with 5.2.

Just curious, Have you tried it on 5.0 ..?? I shall get 5.2 and try running the code on 5.2 and let u knw the details.

Thanks

 

[vatant]

Error

I have a question regarding Mathematica 5.2

I got an error:

NDSolve::ndode: Input is not an ordinary differential equation. More..

It gave me a beep and said, press to "dont show it again" and pressing that gave a result.

If I said "ok", it gave me an serious errror ...

I don't understand if this result is ok ? Since, the error reported by mathematica is "its not an ordinary differential equation" which is true..its a partial differential eqn!

But, I did not get this error when your previous eqn was inputted (with a - sgn before the second term). Could you help explain this?


Thanks !

 

[saltydog]

Hey Vatant,

Can you cut-and-paste your Mathematica code here, then I'll cut-and-paste it back into my 5.2 and run it just in case there is some small syntax error not directly evident. Remember you can't cut-and-paste any special characters so switch the pi's to Pi and other math symbols to character text before the transfer.

 

[vatant]

Hey Vatant,

Can you cut-and-paste your Mathematica code here, then I'll cut-and-paste it back into my 5.2 and run it just in case there is some small syntax error not directly evident. Remember you can't cut-and-paste any special characters so switch the pi's to Pi and other math symbols to character text before the transfer.

Hey there:
The code is essentially the same code you put in. I just copied and pasted it in 5.2 or 5.0 ..and it gives me an error saying the input is not a differential equation and beeps.

Its strange that I dint the error before the change in sign.

In my version of 5.0 or 5.2 your original code with a "-" sign ( for the II term works) while the replacing the "-" with an "+" sign doesn't work.

What do you suggest?

I pressed "dont show the message again" and a plot was obtained similar to yours.

I am confused abt this.

Let me know your thoughts.

Cheers

 

[saltydog]

Hey there:
The code is essentially the same code you put in. I just copied and pasted it in 5.2 or 5.0 ..and it gives me an error saying the input is not a differential equation and beeps.

Its strange that I dint the error before the change in sign.

In my version of 5.0 or 5.2 your original code with a "-" sign ( for the II term works) while the replacing the "-" with an "+" sign doesn't work.

What do you suggest?

I pressed "dont show the message again" and a plot was obtained similar to yours.

I am confused abt this.

Let me know your thoughts.

Cheers

Hey Vatant, the cut-and-paste is not working from PF to Mathematica. Apparently some spurious data/characters are being transmitted, maybe some non-printing escape characters. I get the same problem when I cut-and-paste from here into Mathematica.

Try inputting the NDSolve code from scratch into 5.2 and I bet a dollar it works.

 

[vatant]

Not sure abt the working

I am not sure if the code with a changed "sign" actually works alrite? If you see your solution, with a + and - sign for the second term, the results look similar..and i tried different Mathematica versions to test the code..the code just doesn't work. It says, ND solve not a diff eqn.

I jus copied your code and haven't started working on my own yet. I was thinking of implementing my code on the same guidelines as you suggested, but ran into lots of problems.

Let me know your thoughts.

Can you try implementing some other BC and see if they still work?

Thanks

 

[saltydog]

I am not sure if the code with a changed "sign" actually works alrite? If you see your solution, with a + and - sign for the second term, the results look similar..and i tried different Mathematica versions to test the code..the code just doesn't work. It says, ND solve not a diff eqn.

I jus copied your code and haven't started working on my own yet. I was thinking of implementing my code on the same guidelines as you suggested, but ran into lots of problems.

Let me know your thoughts.

Can you try implementing some other BC and see if they still work?

Thanks

I did try the code with the positive sign and the results seem very similar. However I did back-substitute the results and it did seem to satisfy the differential equation. These are my thoughts:

1. Should have mentioned from the get-go that I'm not at all confident I understand what constitutes a well-posed problem for a forth-order PDE with mixed partials.

2. Perhaps the forth-order term dominates the behavior thus giving rise to the similar results with changing the sign of the second term. When I change the sign of the 4-order term Mathematica buckles.

3. When I change the boundary conditions I too obtain Mathematica errors.

4. I think we're approaching this problem wrong. If it were mine I would go to the library and obtain some book on numerical methods for PDEs preferably with examples of high-ordered ones. First get just one 4-order PDE, anyone preferably with mixed partials working that you know is correct. Then go back to yours.

5. Also I'd seek help in a math department at a chosen university (just go and starting asking questions, you'll get led to the right person). Also, there are reference books on PDEs (look-up tables). I'd check those too. Maybe your equation is addressed there.

 

[vatant]

I did try the code with the positive sign and the results seem very similar. However I did back-substitute the results and it did seem to satisfy the differential equation. These are my thoughts:

1. Should have mentioned from the get-go that I'm not at all confident I understand what constitutes a well-posed problem for a forth-order PDE with mixed partials.

2. Perhaps the forth-order term dominates the behavior thus giving rise to the similar results with changing the sign of the second term. When I change the sign of the 4-order term Mathematica buckles.

3. When I change the boundary conditions I too obtain Mathematica errors.

4. I think we're approaching this problem wrong. If it were mine I would go to the library and obtain some book on numerical methods for PDEs preferably with examples of high-ordered ones. First get just one 4-order PDE, anyone preferably with mixed partials working that you know is correct. Then go back to yours.

5. Also I'd seek help in a math department at a chosen university (just go and starting asking questions, you'll get led to the right person). Also, there are reference books on PDEs (look-up tables). I'd check those too. Maybe your equation is addressed there.

my problem as posted earlier is similar to your example with a "+" sign in front of the second term with the following conditons.

say:

u(0,t)=0
dudt(0,t)=0
dudz(0,t)=-pi
du2dz2(0,t)=0

If atleast your code is working, do you think these boundary conditons give any approrpriat result.

Note that all these values are given at the domain inlet (z,t) where z=0 for any time and evolving in time.

 

[saltydog]

my problem as posted earlier is similar to your example with a "+" sign in front of the second term with the following conditons.

say:

u(0,t)=0
dudt(0,t)=0
dudz(0,t)=-pi
du2dz2(0,t)=0

If atleast your code is working, do you think these boundary conditons give any approrpriat result.

Note that all these values are given at the domain inlet (z,t) where z=0 for any time and evolving in time.

You see, that's where "well-posed" comes in. Your conditions above are only at the left boundary of the domain. What about the right side and what about the initial conditions? As I said above, I'm really not clear on what is sufficient. However I think I know a route to follow: First obtain a good text on PDEs dealing with high ordered partials and IBVPs thereof and see what conditions are used for their solution. What about asking around in a math department for assistance? You have that option?

 

[abra77]

Hi, does anybody know how maple transform the second order time to first order?
Will it be the same as in MATLAB by introducing new variable to reduce the order. Regards