Message Boards Message Boards

0
|
2666 Views
|
2 Replies
|
0 Total Likes
View groups...
Share
Share this post:

Num. Solution of Partial Differential Eq.--The Numerical Method of Lines

Posted 11 years ago
Hi All,

I am trying to solve partial diff. eq. by using method of Numerical Method of Lines

This is the link

https://reference.wolfram.com/mathematica/tutorial/NDSolvePDE.html

I've got some idea

But it does not give me the right graph.

Can anyone check my Mathematica code?

This is my PDE
Subscript[h, t] = 0=
Subscript[\[Alpha], 1] Subscript[h, xxxx] +
  Subscript[\[Alpha], 2] Subscript[h, xx] +
  Subscript[\[Alpha], 5] Subscript[h, xx] Subscript[h, x]^2

I really appriciate it

 n = 100;(*the number of steps*)
 L = 2*\[Pi];(*the length of the period*)
 h = L/n;(*the size of the steps*)
 k = 1;(*frequency*)
 \[Delta] = 0.5
 U[t_] = Table[Subscript[u, i][t], {i, 0, n}];
 initc = Thread[U[0] == Table[0, {n + 1}]];
 Subscript[\[Alpha], 1] = -8; Subscript[\[Alpha], 2] = -216; \
 Subscript[\[Alpha], 5] = -3 Subscript[\[Alpha], 2]/2;
wavenum =
  Sqrt[(Subscript[\[Alpha], 2]/Subscript[\[Alpha], 1])]/Sqrt[2];

count1 = 1;
While[count1 <= n + 1,
  initc[[count1]] =
   Thread[U[0][[count1]] == 1 + \[Delta] Cos[k (count1 - 1) h]];
  count1++];
count1 = 1;
N[initc]
count2 = 1;
eqnsnes = Table[0, {n + 1}];

While[count2 <= n + 1,
  eqnsnes[[count2]] = Thread[D[U[t][[count2]], t] ==
    
     If[count2 == 1,
      Subscript[\[Alpha],
        1] (1/h^4 U[t][[count2 + 2]] - 4 U[t][[count2 + 1]] +
          6 U[t][[count2]] - 4 U[t][[n]] + U[t][[n - 1]]) +
       Subscript[\[Alpha],
        2] ((U[t][[count2 + 1]] - 2 U[t][[count2]] + U[t][[n]])/h^2)/
         wavenum^2 +
       Subscript[\[Alpha],
        5] ((U[t][[count2 + 1]] - 2 U[t][[count2]] + U[t][[n]])/h^2)
        *((U[t][[count2 + 1]] - U[t][[n]])/(2 h))^2, 
     
      If[count2 == 2,
       Subscript[\[Alpha],
         1] (1/h^4 (U[t][[count2 + 2]] - 4 U[t][[count2 + 1]] +
             6 U[t][[count2]] - 4 U[t][[count2 - 1]] + U[t][[n]]) ) +
        Subscript[\[Alpha],
         2] ((U[t][[count2 + 1]] - 2 U[t][[count2]] +
             U[t][[count2 - 1]])/h^2)/wavenum^2 +
       
        Subscript[\[Alpha],
         5] ((U[t][[count2 + 1]] - 2 U[t][[count2]] +
            U[t][[count2 - 1]])/
          h^2)*((U[t][[count2 + 1]] - U[t][[count2 - 1]])/(2 h))^2,
      
       If[count2 == n,
        Subscript[\[Alpha],
          1] (1/h^4 ((U[t][[2]] - 4 U[t][[count2 + 1]] +
               6 U[t][[count2]] - 4 U[t][[count2 - 1]] +
               U[t][[count2 - 2]]) )) +
         Subscript[\[Alpha],
          2] ((U[t][[count2 + 1]] - 2 U[t][[count2]] +
              U[t][[count2 - 1]])/h^2)/wavenum^2 +
         Subscript[\[Alpha],
          5] ((U[t][[count2 + 1]] - 2 U[t][[count2]] +
             U[t][[count2 - 1]])/
           h^2)*((U[t][[count2 + 1]] - U[t][[count2 - 1]])/(2 h))^2,
       
        If[count2 == n + 1,
         Subscript[\[Alpha],
           1] (1/h^4 (U[t][[3]] - 4 U[t][[2]] + 6 U[t][[count2]] -
               4 U[t][[count2 - 1]] + U[t][[count2 - 2]]) ) +
          Subscript[\[Alpha],
           2] ((U[t][[2]] - 2 U[t][[count2]] + U[t][[count2 - 1]])/
             h^2)/wavenum^2 +
          Subscript[\[Alpha],
           5] ((U[t][[2]] - 2 U[t][[count2]] + U[t][[count2 - 1]])/
            h^2)*((U[t][[2]] - U[t][[count2 - 1]])/(2 h))^2,
         Subscript[\[Alpha],
           1] (1/h^4 (U[t][[count2 + 2]] - 4 U[t][[count2 + 1]] +
               6 U[t][[count2]] - 4 U[t][[count2 - 1]] +
               U[t][[count2 - 2]]) ) +
          Subscript[\[Alpha],
           2] ((U[t][[count2 + 1]] - 2 U[t][[count2]] +
               U[t][[count2 - 1]])/h^2)/wavenum^2 +
          Subscript[\[Alpha],
           5] ((U[t][[count2 + 1]] - 2 U[t][[count2]] +
              U[t][[count2 - 1]])/
            h^2)*((U[t][[count2 + 1]] - U[t][[count2 - 1]])/(
            2 h))^2]]]]]; count2++];
eqnsnes

lines = NDSolve[{eqnsnes, initc}, U[t], {t, 0, 1}, MaxSteps -> 20000]

ParametricPlot3D[
Evaluate[Table[{i*h, t, First[Subscript[u, i][t] /. lines]}, {i, 0,
    n}]], {t, 0, 1}, PlotRange -> All,
AxesLabel -> {"Spatial(x)", "Time(t)", "Height(h)"}]
POSTED BY: selahittin cinar
2 Replies
yes, it is second derivative. I dont why it shows that way. Also, in my code variable is t instead of x..

how can I copy and paste my code as it is orijinal?
POSTED BY: selahittin cinar
Are you thinking that Subscript[h, xx] performs a differentiation? 
It is just a static typeset thing.
POSTED BY: Bruce Miller
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract