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)"}]