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

1. Jan 27, 2006

### 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

Last edited: Jan 27, 2006
2. Jan 27, 2006

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

3. Jan 28, 2006

### 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:

Code (Text):

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:

#### Attached Files:

• ###### wave eq.JPG
File size:
18.8 KB
Views:
104
Last edited: Jan 28, 2006
4. Jan 28, 2006

### 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:

Code (Text):

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

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:

Code (Text):

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.

#### Attached Files:

• ###### dif of wave solution.JPG
File size:
5.3 KB
Views:
69
Last edited: Jan 28, 2006
5. Jan 29, 2006

### 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:

Code (Text):

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.

File size:
18.5 KB
Views:
82
File size:
10.4 KB
Views:
68
6. Feb 2, 2006

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

7. Feb 2, 2006

### saltydog

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

Code (Text):

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.

8. Feb 2, 2006

### 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 doesnt acknowlege that !! did u try that sign change? says not a Diff eqn..what do you suggest?

9. Feb 3, 2006

### saltydog

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?

10. Feb 3, 2006

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

Last edited: Feb 3, 2006
11. Feb 3, 2006

### saltydog

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

Code (Text):

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?

#### Attached Files:

• ###### +pde.JPG
File size:
19 KB
Views:
72
12. Feb 3, 2006

### vatant

Mathematica 5.0

I am using version 5.0 of mathematica and looks like it doesnt gonna work..you think thats 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

13. Feb 4, 2006

### 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 dont 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 !

14. Feb 4, 2006

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

15. Feb 4, 2006

### vatant

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

Cheers

16. Feb 5, 2006

### saltydog

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.

17. Feb 6, 2006

### 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 doesnt work. It says, ND solve not a diff eqn.

I jus copied your code and havent 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.

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

Thanks

18. Feb 7, 2006

### saltydog

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, any one 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.

19. Feb 8, 2006

### vatant

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.

20. Feb 9, 2006

### saltydog

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?