0
$\begingroup$

I want to plot the solution obtained from NDSolve as a multicurve plot like shown in the following sample image (image ref). Here, x1, x2 , x3, x4 .. are the solution computed at includePoints = {{10}, {20}, {30}, {40}, {50}} in my case.

enter image description here

In the following code (ref):

Needs["NDSolve`FEM`"]
region = Line[{{0}, {100}}];
includePoints = {{10}, {20}, {30}, {40}, {50}};
mesh = ToElementMesh[region, "IncludePoints" -> includePoints, 
  "MaxCellMeasure" -> 0.0008]
vars = {c[t, x], t, {x}};
RegularizedDeltaPoint[g_, X_List, Xs_List] := 
 Piecewise[{{Times @@ Thread[1/(4 g) (1 + Cos[\[Pi]/(2 g) (X - Xs)])],
     And @@ Thread[RealAbs[X - Xs] <= 2 g]}, {0, True}}]
Subscript[h, mesh] = Sqrt[Min[mesh["MeshElementMeasure"]]];
Subscript[gamma, reg] = Subscript[h, mesh]/2;
temp = RegularizedDeltaPoint[Subscript[gamma, reg], {x}, 
   includePoints[[1]]];
parameters = {kappa -> {{910}}, v1 -> 162, 
   gamma -> Subscript[gamma, reg], Qp -> 1.5};
pde = {Derivative[1, 0][c][t, x] + 
      Inactive[
        Div][(-kappa).Inactive[Grad][
         c[t, x], {x}], {x}] + {v1}.Inactive[Grad][c[t, x], {x}] + 
      Qp*RegularizedDeltaPoint[gamma, {x}, {10}] == 0, 
    c[0, x] == 1} /. parameters;

tEnd = 2;
cfun = NDSolveValue[{pde, DirichletCondition[c[t, x] == 5, x == 0]}, c, {t, 0, tEnd}, {x} \[Element] mesh];
Manipulate[
 Plot[cfun[t, x], {x} \[Element] region, 
  PlotRange -> {{0, 100}, {0, 5}}], {t, 0, tEnd}]

enter image description here

The above figure displays solution vs x-position. Instead, I want to plot solution curves observed at x-positions in includePoints = {{10}, {20}, {30}, {40}, {50}}; as a function of time. as a function of time.

Suggestions will be really appreciated.

EDIT:

  1. May I also know how to save the solution in cfun at includePoints = {{10}, {20}, {30}, {40}, {50}} over the integration time span to a text file?

  2. Using the following command

    With[{i = Flatten[includePoints]}, Plot[Evaluate[cfun[t, #] & /@ i], {t, 0, tEnd}, PlotRange -> {{0, tEnd}, {0, 5}}, PlotLegends -> (StringTemplate["c(t,``)"] /@ i)]]

I could generate, enter image description here

But I am not sure why c(t,0), the first point which is the left boundary isn't at 5 (Dirichlet conditioned defined) at t=0. Could someone please have a look? Also, the legend labels aren't formatted correctly String Template appears as text in the legend).

$\endgroup$
4
  • $\begingroup$ @Michael E2 Could you please help me with this? $\endgroup$
    – Natasha
    Jun 20 at 7:18
  • $\begingroup$ Try including the Dirichlet condition in a list with the PDE, like this cfun = NDSolveValue[{pde, DirichletCondition[c[t, x] == 5, x == 0]}, c, {t, 0, tEnd}, {x} \[Element] mesh]; $\endgroup$
    – LouisB
    Jun 20 at 8:11
  • $\begingroup$ @LouisB Thank you, I tried the above but it is still not clear to me how the plot can be generated. Could you please explain a bit more? Or, may I know how to save the solution in cfun at includePoints = {{10}, {20}, {30}, {40}, {50}} over the integration time span to a text file? $\endgroup$
    – Natasha
    Jun 20 at 10:40
  • $\begingroup$ This question is a continuation of a question here $\endgroup$
    – user21
    Jun 21 at 4:31
5
$\begingroup$

Update to fix boundary and initial condition inconsistency

Your "MaxCellMeasure" is much more refined than it needs to be. Furthermore, you will note that the initial condition and the DirichletCondition are inconsistent. This causes the DirichletCondition to be ignored. One way to remove this inconsistency is to rapidly ramp up the DirichletCondition from the initial condition to its desired value.

Needs["NDSolve`FEM`"]
region = Line[{{0}, {100}}];
includePoints = {{10}, {20}, {30}, {40}, {50}};
mesh = ToElementMesh[region, "IncludePoints" -> includePoints, 
   "MaxCellMeasure" -> 1];
vars = {c[t, x], t, {x}};
RegularizedDeltaPoint[g_, X_List, Xs_List] := 
 Piecewise[{{Times @@ Thread[1/(4 g) (1 + Cos[π/(2 g) (X - Xs)])],
     And @@ Thread[RealAbs[X - Xs] <= 2 g]}, {0, True}}]
Subscript[h, mesh] = Sqrt[Min[mesh["MeshElementMeasure"]]];
Subscript[gamma, reg] = Subscript[h, mesh]/2;
temp = RegularizedDeltaPoint[Subscript[gamma, reg], {x}, 
   includePoints[[1]]];
parameters = {kappa -> {{910}}, v1 -> 162, 
   gamma -> Subscript[gamma, reg], Qp -> 1.5};
pde = {Derivative[1, 0][c][t, x] + 
      Inactive[Div][(-kappa) . 
        Inactive[Grad][c[t, x], {x}], {x}] + {v1} . 
       Inactive[Grad][c[t, x], {x}] + 
      Qp*RegularizedDeltaPoint[gamma, {x}, {10}] == 0, c[0, x] == 1} /. 
   parameters;

tEnd = 2;
cfun = NDSolveValue[
   pde~Join~{DirichletCondition[c[t, x] == 4 (1 - Exp[-1000 t]) + 1, 
      x == 0]}, c, {t, 0, tEnd}, {x} ∈ mesh];
With[{i = Flatten[{0}~Join~includePoints]}, 
 Plot[Evaluate[cfun[t, #] & /@ i], {t, 0, tEnd}, 
  PlotRange -> {{0, tEnd}, {0, 5.1}}, 
  PlotLegends -> (StringTemplate["c(t,``)"] /@ i)]]
ListPlot[Table[cfun[t, x], {x, 0, 100, 10}, {t, 0, tEnd, tEnd/100}], 
 Joined -> True, DataRange -> {0, tEnd}]

Solution with fixed inconsistency

Original answer

Here is one way you could plot the solutions:

With[{i = Flatten[includePoints]},
 Plot[Evaluate[cfun[t, #] & /@ i], {t, 0, tEnd}, 
  PlotRange -> {{0, tEnd}, {0, 5}}, 
  PlotLegends -> (StringTemplate["c(t,``)"] /@ i)]
 ]

enter image description here

$\endgroup$
5
  • $\begingroup$ Thanks so much, works great. Sorry for another question. Could you please have a look at my edit? $\endgroup$
    – Natasha
    Jun 20 at 12:44
  • $\begingroup$ @Natasha If I quit the kernel and copy the code from your post, I get the image shown in my answer. Unless you changed the definition of includePoints, 0 should not be there. What $Version of Mathematica are you using? $\endgroup$
    – Tim Laska
    Jun 20 at 13:04
  • $\begingroup$ Thanks a ton for the answer :) . I'm using version 12.0 $\endgroup$
    – Natasha
    Jun 20 at 17:11
  • $\begingroup$ Could you please have a look at this post? $\endgroup$
    – Natasha
    Jun 21 at 16:18
  • $\begingroup$ Is there an option to specify MinCellMeasure ? Could you please have a look at this post? $\endgroup$
    – Natasha
    Jun 28 at 1:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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