Group Abstract

Message Boards

WOLFRAM COMMUNITY

6.5K Views

16 Replies

0 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Solve the following delay differential equations?

Rachele Allena

Posted 7 years ago

Dear all, I am solving a PDE such as : Button["Stop", stop = True] stop = False; currentTime = "Initialization"; interval = T; SetOptions[EvaluationNotebook[], WindowStatusArea -> Dynamic["t = " <> ToString[CForm[currentTime]]]]; {uif, vif} = NDSolveValue[{Activate[divsig - divsig0 == inertia + fadh], ci, cid, WhenEvent[stop, end = t; "StopIntegration"]}, {u, v}, {x, y} \[Element] mesh, {t, 0, interval}, Method -> { "PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MeshOrder" -> 2, "MaxCellMeasure" -> 0.5}, "IntegrationOrder" -> 4}}}, EvaluationMonitor :> (currentTime = t)(* ,AccuracyGoal\[Rule]5*), MaxStepSize -> 0.5]; Now, the variable 'divsig0' depends on a previous defined variable, 'eaf', which in turn depends itself on the solution of the PDE through 'uif' like this eaf = Table[ If[pcontour[[i, j]][[1]] + uif[t, pcontour[[i, j]][[1]], pcontour[[i, j]][[2]]] > pillarscenters[[i, j]][[1]], t/100 , 0], {i, Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}]; It seems I should used a Delay Differential Equation according to the documentation, but I have trouble to specify the initial hystory function...could anyone help me? I can send the whole code if necessary... Best RA

Dear all,

I am solving a PDE such as :

Button["Stop", stop = True]
stop = False;
currentTime = "Initialization";
interval = T;
SetOptions[EvaluationNotebook[], 
  WindowStatusArea -> Dynamic["t = " <> ToString[CForm[currentTime]]]];
{uif, vif} = 
  NDSolveValue[{Activate[divsig - divsig0 == inertia + fadh], ci, cid,
     WhenEvent[stop, end = t; "StopIntegration"]}, {u, 
    v}, {x, y} \[Element] mesh, {t, 0, interval}, Method -> {
     "PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement", 
         "MeshOptions" -> {"MeshOrder" -> 2, "MaxCellMeasure" -> 0.5},
          "IntegrationOrder" -> 4}}}, 
   EvaluationMonitor :> (currentTime = t)(* ,AccuracyGoal\[Rule]5*), 
   MaxStepSize -> 0.5];

Now, the variable 'divsig0' depends on a previous defined variable, 'eaf', which in turn depends itself on the solution of the PDE through 'uif' like this

eaf = Table[
   If[pcontour[[i, j]][[1]] + 
     uif[t, pcontour[[i, j]][[1]], pcontour[[i, j]][[2]]] > pillarscenters[[i, j]][[1]], t/100 , 0], {i, 
    Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}];

It seems I should used a Delay Differential Equation according to the documentation, but I have trouble to specify the initial hystory function...could anyone help me?

I can send the whole code if necessary...

Best RA

POSTED BY: Rachele Allena

16 Replies

Sort By:

Rachele Allena

Posted 7 years ago

'eaf' is an initial deformation which provide 'divsig0' According to its expression eaf[kk_] := Table[If[pcontour[[i]][[1]] + uif[kk][t, pcontour[[i]][[1]], pcontour[[i]][[2]]] > pillarscenters[[i]][[1]], t/100, 0], {i, Dimensions[alpha][[1]]}]; I would like that such an initial deformation stops once the the sum pcontour[[i]][[1]] + uif[kk][t, pcontour[[i]][[1]], pcontour[[i]][[2]]] is less than pillarscenters[[i]][[1]] By doing so, when plotting the displacement of a point in the mesh, this should increase until a certain time point and than stay constant, whereas when plotting 'eaf' this should increase until the same time point and than drop to zero according to the If function However, when I do 1 iteration it works for eaf but not for the displacement which continously increase over time and I do not understand why since if I test the condition pcontour[[i]][[1]] + uif[kk][t, pcontour[[i]][[1]], pcontour[[i]][[2]]] > pillarscenters[[i]][[1]] it evaluates to true for the first time steps and false for the rest of the interval, than the displacement should be stop too...

POSTED BY: Rachele Allena

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 7 years ago

What problem are you solving? What do all these `If[]` mean?

POSTED BY: Alexander Trounev

Rachele Allena

Posted 7 years ago

Actually it runs like this, but I am not sure that this war 'uif' is correctly taken into accoutn then. Normally, eaf[kk] should be equal to 0, and therefore the displacement, once pcontour[[i, j]][[1]] + uif[kk][t, pcontour[[i, j]][[1]], pcontour[[i, j]][[2]]] > pillarscenters[[i, j]][[1]] is not satisfied anymore. But when you plot Plot[{uif[kf][t, pcontour[[1, 1]][[1]], pcontour[[1, 1]][[2]]], vif[kf][t, pcontour[[1, 1]][[1]], pcontour[[1, 1]][[2]]]}, {t, 0, tf}, PlotRange -> All] this does not occur. Then I am wondering if the 'If' in the definition of 'eaf[kk] is correctly taken into account? best rachele

POSTED BY: Rachele Allena

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 7 years ago

No, you're not right, we made this replacement in `(divsig0 /. uif -> uif[i - 1])` in line `Do[]`. If you want to make this replacement earlier, then make it so that there is no error due to the coincidence of the names of the iterator eaf[kk_] := Table[If[pcontour[[i, j]][[1]] + uif[kk][t, pcontour[[i, j]][[1]], pcontour[[i, j]][[2]]] > pillarscenters[[i, j]][[1]], t/100, 0], {i, Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}]; ear = Table[{eax[[i, j]] ea c1hs[ea, 10^(-6)], eay[[i, j]] ea c1hs[ea, 10^(-6)]}, {i, Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}]; mesh = ToElementMesh[cell, MaxCellMeasure -> {"Length" -> 1 10^(-6)}, "BoundaryMeshGenerator" -> "Continuation"]; Ecell = 10^4; EY1 = Ecell(*rigid+(Ecell/10)soft); EY2 = Ecell/7; nu12 = 0.3; nu21 = 0.3; G12 = EY1/(2(1 + nu12)); ro = 1000; lzonef = 210^(-6); lzoner = -210^(-6); zonef = c1hs[(x nppodx[[1, 1]] + y nppody[[1, 1]]) - lzonef, 10^(-12)]; zoner = c1hs[-(x nppodx[[1, 1]] + y nppody[[1, 1]]) + lzoner, 10^(-12)]; muv = 10^(16); mustab = 10^(3); mupil = 10^(9); fadhf = muvzonef c1hs[-dea, 10^(-10)] {D[u[t, x, y], t], D[v[t, x, y], t]}; fadhr = muvzoner c1hs[dea, 10^(-10)] {D[u[t, x, y], t], D[v[t, x, y], t]}; fstab = mustab {D[u[t, x, y], t], D[v[t, x, y], t]}; fpillars = mupil {D[u[t, x, y], t], D[v[t, x, y], t]}soft; fadh = fadhf + fadhr + fstab; mechprop = {Y -> EY1, \[Nu] -> nu12, rho -> ro}; \[Lambda]dp = Y \[Nu]/(1 + \[Nu])/(1 - 2 \[Nu]); \[Lambda] = Y \[Nu]/(1 - \[Nu]^2); \[Mu] = Y/2/(1 + \[Nu]); \[CurlyEpsilon] = {{ix.Inactive[Grad][u[t, x, y], {x, y}], 1/2 (iy.Inactive[Grad][u[t, x, y], {x, y}] + ix.Inactive[Grad][v[t, x, y], {x, y}])}, {1/ 2 (iy.Inactive[Grad][u[t, x, y], {x, y}] + ix.Inactive[Grad][v[t, x, y], {x, y}]), iy.Inactive[Grad][v[t, x, y], {x, y}]}}; \[CurlyEpsilon]0[ kk_] := -ppods[[1, 1]] (eaf[kk][[1, 1]]) ((x nppodx[[1, 1]] + y nppody[[1, 1]])/ 10^(-6))^2 {{nppodx[[1, 1]] nppodx[[1, 1]], nppodx[[1, 1]] nppody[[1, 1]]}, {nppodx[[1, 1]] nppody[[1, 1]], nppody[[1, 1]] nppody[[1, 1]]}}; \[Sigma] = \[Lambda] Tr[\[CurlyEpsilon]] Id + 2 \[Mu] (\[CurlyEpsilon]); \[Sigma]0[kk_] := \[Lambda] Tr[\[CurlyEpsilon]0[kk]] Id + 2 \[Mu] (\[CurlyEpsilon]0[kk]); divsig = {Inactive[Div][(\[Sigma].ix), {x, y}], Inactive[Div][(\[Sigma].iy), {x, y}]} /. mechprop; divsig0[kk_] := {Inactive[Div][(\[Sigma]0[kk].ix), {x, y}], Inactive[Div][(\[Sigma]0[kk].iy), {x, y}]} /. mechprop; inertia = rho {D[u[t, x, y], {t, 2}], D[v[t, x, y], {t, 2}]} /. mechprop; ci = {u[0, x, y] == 0., v[0, x, y] == 0.}; cid = {Derivative[1, 0, 0][u][0, x, y] == 0, Derivative[1, 0, 0][v][0, x, y] == 0}; uif[0][t_, x_, y_] := 0 vif[0][t_, x_, y_] := 0 kf = 4; tf = 15; Do[{uif[kk], vif[kk]} = NDSolveValue[{Activate[divsig - divsig0[kk - 1] == inertia + fadh], ci, cid}, {u, v}, {x, y} \[Element] mesh, {t, 0, tf}, Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MeshOrder" -> 2, "MaxCellMeasure" -> 0.5}, "IntegrationOrder" -> 4}}}, MaxStepSize -> 0.5];, {kk, 1, kf}] Table[ContourPlot[uif[i][tf, x, y], {x, y} \[Element] mesh, PlotLegends -> Automatic, Contours -> 20, ColorFunction -> "TemperatureMap", PlotLabel -> Row[{"i=", i}]], {i, 1, kf}] Attachments:* uv4N.nb

No, you're not right, we made this replacement in (divsig0 /. uif -> uif[i - 1]) in line Do[]. If you want to make this replacement earlier, then make it so that there is no error due to the coincidence of the names of the iterator

eaf[kk_] := 
  Table[If[pcontour[[i, j]][[1]] + 
      uif[kk][t, pcontour[[i, j]][[1]], pcontour[[i, j]][[2]]] > 
     pillarscenters[[i, j]][[1]], t/100, 0], {i, 
    Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}];
ear = Table[{eax[[i, j]] ea c1hs[ea, 10^(-6)], 
    eay[[i, j]] ea c1hs[ea, 10^(-6)]}, {i, 
    Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}];

mesh = ToElementMesh[cell, MaxCellMeasure -> {"Length" -> 1 10^(-6)}, 
   "BoundaryMeshGenerator" -> "Continuation"];

Ecell = 10^4;
EY1 = Ecell(**rigid+(Ecell/10)*soft*);
EY2 = Ecell/7;
nu12 = 0.3;
nu21 = 0.3;
G12 = EY1/(2*(1 + nu12));
ro = 1000;

lzonef = 2*10^(-6);
lzoner = -2*10^(-6);
zonef = c1hs[(x nppodx[[1, 1]] + y nppody[[1, 1]]) - lzonef, 10^(-12)];
zoner = c1hs[-(x nppodx[[1, 1]] + y nppody[[1, 1]]) + lzoner, 
   10^(-12)];
muv = 10^(16);
mustab = 10^(3);
mupil = 10^(9);
fadhf = muv*zonef*
   c1hs[-dea, 10^(-10)] {D[u[t, x, y], t], D[v[t, x, y], t]};
fadhr = muv*zoner*
   c1hs[dea, 10^(-10)] {D[u[t, x, y], t], D[v[t, x, y], t]};
fstab = mustab {D[u[t, x, y], t], D[v[t, x, y], t]};
fpillars = mupil {D[u[t, x, y], t], D[v[t, x, y], t]}*soft;
fadh = fadhf + fadhr + fstab;

mechprop = {Y -> EY1, \[Nu] -> nu12, rho -> ro};
\[Lambda]dp = Y \[Nu]/(1 + \[Nu])/(1 - 2 \[Nu]);
\[Lambda] = Y \[Nu]/(1 - \[Nu]^2);
\[Mu] = Y/2/(1 + \[Nu]);
\[CurlyEpsilon] = {{ix.Inactive[Grad][u[t, x, y], {x, y}], 
    1/2 (iy.Inactive[Grad][u[t, x, y], {x, y}] + 
       ix.Inactive[Grad][v[t, x, y], {x, y}])}, {1/
      2 (iy.Inactive[Grad][u[t, x, y], {x, y}] + 
       ix.Inactive[Grad][v[t, x, y], {x, y}]), 
    iy.Inactive[Grad][v[t, x, y], {x, y}]}};
\[CurlyEpsilon]0[
   kk_] := -ppods[[1, 
     1]] (eaf[kk][[1, 1]]) ((x nppodx[[1, 1]] + y nppody[[1, 1]])/
      10^(-6))^2 {{nppodx[[1, 1]] nppodx[[1, 1]], 
     nppodx[[1, 1]] nppody[[1, 1]]}, {nppodx[[1, 1]] nppody[[1, 1]], 
     nppody[[1, 1]] nppody[[1, 1]]}};

\[Sigma] = \[Lambda] Tr[\[CurlyEpsilon]] Id + 
   2 \[Mu] (\[CurlyEpsilon]);
\[Sigma]0[kk_] := \[Lambda] Tr[\[CurlyEpsilon]0[kk]] Id + 
   2 \[Mu] (\[CurlyEpsilon]0[kk]);
divsig = {Inactive[Div][(\[Sigma].ix), {x, y}], 
    Inactive[Div][(\[Sigma].iy), {x, y}]} /. mechprop;
divsig0[kk_] := {Inactive[Div][(\[Sigma]0[kk].ix), {x, y}], 
    Inactive[Div][(\[Sigma]0[kk].iy), {x, y}]} /. mechprop;
inertia = 
  rho {D[u[t, x, y], {t, 2}], D[v[t, x, y], {t, 2}]} /. mechprop;
ci = {u[0, x, y] == 0., v[0, x, y] == 0.};
cid = {Derivative[1, 0, 0][u][0, x, y] == 0, 
   Derivative[1, 0, 0][v][0, x, y] == 0};
uif[0][t_, x_, y_] := 0
vif[0][t_, x_, y_] := 0
kf = 4; tf = 15;

Do[{uif[kk], vif[kk]} = 
   NDSolveValue[{Activate[divsig - divsig0[kk - 1] == inertia + fadh],
      ci, cid}, {u, v}, {x, y} \[Element] mesh, {t, 0, tf}, 
    Method -> {"PDEDiscretization" -> {"MethodOfLines", 
        "SpatialDiscretization" -> {"FiniteElement", 
          "MeshOptions" -> {"MeshOrder" -> 2, 
            "MaxCellMeasure" -> 0.5}, "IntegrationOrder" -> 4}}}, 
    MaxStepSize -> 0.5];, {kk, 1, kf}]
Table[ContourPlot[uif[i][tf, x, y], {x, y} \[Element] mesh, 
  PlotLegends -> Automatic, Contours -> 20, 
  ColorFunction -> "TemperatureMap", PlotLabel -> Row[{"i=", i}]], {i,
   1, kf}]

fig3

POSTED BY: Alexander Trounev

Rachele Allena

Posted 7 years ago

Actually it works like this, but the problem is that the line eaf = Table[ If[pcontour[[i, j]][[1]] + uif[t, pcontour[[i, j]][[1]], pcontour[[i, j]][[2]]] > pillarscenters[[i, j]][[1]], t/100, 0], {i, Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}]; is wrong because one needs to write eaf = Table[ If[pcontour[[i, j]][[1]] + uif[[kf]t, pcontour[[i, j]][[1]], pcontour[[i, j]][[2]]] > pillarscenters[[i, j]][[1]], t/100, 0], {i, Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}]; since with the iteration method, uif depends now on kf too. Am I right? Then withis change, it does not work anymore...

Actually it works like this, but the problem is that the line

eaf = Table[
   If[pcontour[[i, j]][[1]] + 
      uif[t, pcontour[[i, j]][[1]], pcontour[[i, j]][[2]]] > 
     pillarscenters[[i, j]][[1]], t/100, 0], {i, 
    Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}];

is wrong because one needs to write

 eaf = Table[
       If[pcontour[[i, j]][[1]] + 
          uif[[kf]t, pcontour[[i, j]][[1]], pcontour[[i, j]][[2]]] > 
         pillarscenters[[i, j]][[1]], t/100, 0], {i, 
        Dimensions[alpha][[1]]}, {j, Dimensions[alpha][[2]]}];

since with the iteration method, uif depends now on kf too. Am I right? Then withis change, it does not work anymore...

POSTED BY: Rachele Allena

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 7 years ago

Check `$Version`. My version "11.3.0 for Microsoft Windows (64-bit) (March 7, 2018)" works good. I calculated 8 iterations. The process branches at `t>10.4`. I attach the file. Attachments: uv8.nb

POSTED BY: Alexander Trounev

Rachele Allena

Posted 7 years ago

POSTED BY: Rachele Allena

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 7 years ago

Remove all unnecessary. The last lines of code should be such ci = {u[0, x, y] == 0., v[0, x, y] == 0.}; cid = {Derivative[1, 0, 0][u][0, x, y] == 0, Derivative[1, 0, 0][v][0, x, y] == 0}; uif[0][t_, x_, y_] := 0 vif[0][t_, x_, y_] := 0 kf = 2; tf = 15; Do[{uif[i], vif[i]} = NDSolveValue[{Activate[ divsig - (divsig0 /. uif -> uif[i - 1]) == inertia + fadh], ci, cid}, {u, v}, {x, y} \[Element] mesh, {t, 0, tf}, Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MeshOrder" -> 2, "MaxCellMeasure" -> 0.5}, "IntegrationOrder" -> 4}}}, MaxStepSize -> 0.5];, {i, 1, kf}] Run the code and make sure the solution on the first and second iteration is different at `t>10.4`. Therefore it is necessary to increase `kf`. Table[ContourPlot[uif[kf][t, x, y], {x, y} \[Element] mesh, PlotLegends -> Automatic, ColorFunction -> "TemperatureMap", Contours -> 20, PlotLabel -> Row[{"t=", t}]], {t, 2, tf, 2}] Plot[Evaluate[Table[uif[i][t, -210^-6, -210^-6], {i, 1, kf}]], {t, 0, tf}, AxesLabel -> {"t", "u"}]

Remove all unnecessary. The last lines of code should be such

ci = {u[0, x, y] == 0., v[0, x, y] == 0.};
cid = {Derivative[1, 0, 0][u][0, x, y] == 0, 
   Derivative[1, 0, 0][v][0, x, y] == 0};
uif[0][t_, x_, y_] := 0
vif[0][t_, x_, y_] := 0
kf = 2; tf = 15;
Do[{uif[i], vif[i]} = 
   NDSolveValue[{Activate[
      divsig - (divsig0 /. uif -> uif[i - 1]) == inertia + fadh], ci, 
     cid}, {u, v}, {x, y} \[Element] mesh, {t, 0, tf}, 
    Method -> {"PDEDiscretization" -> {"MethodOfLines", 
        "SpatialDiscretization" -> {"FiniteElement", 
          "MeshOptions" -> {"MeshOrder" -> 2, 
            "MaxCellMeasure" -> 0.5}, "IntegrationOrder" -> 4}}}, 
    MaxStepSize -> 0.5];, {i, 1, kf}]

Run the code and make sure the solution on the first and second iteration is different at t>10.4. Therefore it is necessary to increase kf.

Table[ContourPlot[uif[kf][t, x, y], {x, y} \[Element] mesh, 
  PlotLegends -> Automatic, ColorFunction -> "TemperatureMap", 
  Contours -> 20, PlotLabel -> Row[{"t=", t}]], {t, 2, tf, 2}]
Plot[Evaluate[Table[uif[i][t, -2*10^-6, -2*10^-6], {i, 1, kf}]], {t, 
  0, tf}, AxesLabel -> {"t", "u"}]

fig2

POSTED BY: Alexander Trounev

Rachele Allena

Posted 7 years ago

POSTED BY: Rachele Allena

Rachele Allena

Posted 7 years ago

This normally shouldn't cause any problem for the NDSolve command, it is just a warning and it is additionally not used in the NDSolve

POSTED BY: Rachele Allena

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 7 years ago

POSTED BY: Alexander Trounev

Rachele Allena

Posted 7 years ago

It is bizarre, when I open the file from scratch, it is not working anymore... I'll try to fix it, but actually even when it worked, it was very slow...was it the same on your machine? Thanks rachele

POSTED BY: Rachele Allena

Rachele Allena

Posted 7 years ago

Hello Thank you very much for the suggestion. However, when I run the code, it computes correctly up to tf = 10 but it does not stop and cannot post process the results. Do you have any idea why? Thanks rachele

POSTED BY: Rachele Allena

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 7 years ago

Here you can use the iteration method. I checked that the solution quickly converged. An example with three iterations uif[0][t_, x_, y_] := 0 vif[0][t_, x_, y_] := 0 kf = 3; tf = 10; Do[{uif[i], vif[i]} = NDSolveValue[{Activate[ divsig - (divsig0 /. uif -> uif[i - 1]) == inertia + fadh], ci, cid}, {u, v}, {x, y} \[Element] mesh, {t, 0, tf}, Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MeshOrder" -> 2, "MaxCellMeasure" -> 0.5}, "IntegrationOrder" -> 4}}}, MaxStepSize -> 0.5];, {i, 1, kf}] Table[ContourPlot[uif[kf][t, x, y], {x, y} \[Element] mesh, PlotLegends -> Automatic, ColorFunction -> "TemperatureMap", Contours -> 20, PlotLabel -> Row[{"t=", t}]], {t, 2, tf, 2}]

Here you can use the iteration method. I checked that the solution quickly converged. An example with three iterations

uif[0][t_, x_, y_] := 0
vif[0][t_, x_, y_] := 0
kf = 3; tf = 10;

Do[{uif[i], vif[i]} = 
   NDSolveValue[{Activate[
      divsig - (divsig0 /. uif -> uif[i - 1]) == inertia + fadh], ci, 
     cid}, {u, v}, {x, y} \[Element] mesh, {t, 0, tf}, 
    Method -> {"PDEDiscretization" -> {"MethodOfLines", 
        "SpatialDiscretization" -> {"FiniteElement", 
          "MeshOptions" -> {"MeshOrder" -> 2, 
            "MaxCellMeasure" -> 0.5}, "IntegrationOrder" -> 4}}}, 
    MaxStepSize -> 0.5];, {i, 1, kf}]

Table[ContourPlot[uif[kf][t, x, y], {x, y} \[Element] mesh, 
  PlotLegends -> Automatic, ColorFunction -> "TemperatureMap", 
  Contours -> 20, PlotLabel -> Row[{"t=", t}]], {t, 2, tf, 2}]

fig1