1
$\begingroup$

Given the analytical solution

f[x_] := Exp[-x^2/2] (1 + 0.2*HermiteH[4, x] + 0.01*HermiteH[6, x])
n = 10;
A = Table[
   1/(Sqrt[Pi] 2^m m!) Integrate[
     f[x] HermiteH[m, x] Exp[x^2/2], {x, 0, 1}], {m, 1, n}];
u[x_, t_] := 
 Sum[A[[m]] Exp[-m/2 t] Exp[-x^2/2] HermiteH[m - 1, x], {m, 1, n}]

and the numerical solution

Clear["Global`*"]

eigen[n_] := n + 1/2
Her[n_] := HermiteH[n, x]
Y[n_] := Exp[-x^2/2]*Her[n]
Y[1]
Y[4]
Y[6]
f[x_] = Y[1] + 0.2 Y[4] + 0.01 Y[6]
IC = u[x, 0] == f[x]
PDE = \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[x, t]\)\) - 1/2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(x, x\)]\(u[x, t]\)\) + 
   1/2*x^2*u[x, t] == 0
BCs = {u[-21, t] == 0, u[21, t] == 0}
sol = First[NDSolve[{PDE, BCs, IC}, u, {t, 0, 10}, {x, -21, 21}]]

How could I compare the two solutions with the plot command? I face some difficulties in this particular case

$\endgroup$
4
  • 1
    $\begingroup$ How you know that u is exact solution for this PDE? $\endgroup$ Mar 28, 2023 at 5:22
  • $\begingroup$ Plot the difference. $\endgroup$ Mar 28, 2023 at 15:14
  • $\begingroup$ @DanielHuber How could I do that, please? $\endgroup$ Mar 28, 2023 at 20:23
  • 1
    $\begingroup$ Look at my answer. $\endgroup$ Mar 29, 2023 at 8:11

1 Answer 1

1
$\begingroup$

Do not use the same symbol for the analytical and numerical solution. Then move the "Clear" to the beginning, otherwise you will clear the numerical solution. Then store a function in "sol" not a rule.

Clear["Global`*"]
f[x_] := Exp[-x^2/2] (1 + 0.2*HermiteH[4, x] + 0.01*HermiteH[6, x])
n = 10;
A = Table[
   1/(Sqrt[Pi] 2^m m!) Integrate[
     f[x] HermiteH[m, x] Exp[x^2/2], {x, 0, 1}], {m, 1, n}];
ua[x_, t_] := 
 Sum[A[[m]] Exp[-m/2 t] Exp[-x^2/2] HermiteH[m - 1, x], {m, 1, n}]

Her[n_] := HermiteH[n, x];
Y[n_] := Exp[-x^2/2]*Her[n];
f[x_] = Y[1] + 0.2 Y[4] + 0.01 Y[6];
IC = u[x, 0] == f[x];
PDE = \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[x, t]\)\) - 1/2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(x, x\)]\(u[x, t]\)\) + 
       1/2*x^2*u[x, t] == 0;
BCs = {u[-21, t] == 0, u[21, t] == 0};
sol[x_, t_] = 
 u[x, t] /. 
  First[NDSolve[{PDE, BCs, IC}, u, {t, 0, 10}, {x, -21, 21}]]

Plot3D[{sol[x, t] - ua[x, t]}, {x, -21, 21}, {t, 0, 10}, 
 PlotRange -> All]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.