## Postprocessing the results of NDSolve for analysis of structures

So, I get the deformations on solving a model of a structure using NDSolve. Now, does Mathematica have post-processors to retrieve the strains and stresses or is it up to me to implement that part?

I could not find anything about retrieving the stresses in the documentation. Am I missing something?

## How to obtain the same result of the same system if I write the system in different ways (NDSolve)?

I did program where I solved a system when I organized in matrix form, this is the code

Clear["Global*"]  SeedRandom[1234]  Nmax = 5; (*Number of sites*)  tini = 0; (*initial time*)  tmax = 200; (*maximal time*)  \[Sigma]2 = 0; (*Variance*)  n0 = 5; (*initial condition*)  ra = 1; (*coupling range*)  \[Psi]ini = Table[KroneckerDelta[n0 - i], {i, 1, Nmax}];  RR = RandomReal[{-Sqrt[3*\[Sigma]2], Sqrt[3*\[Sigma]2]}, Nmax];  Z = Table[     Sum[KroneckerDelta[i - j + k], {k, 1, ra}] +       Sum[KroneckerDelta[i - j - k], {k, 1, ra}], {i, 1, Nmax}, {j, 1,       Nmax}] + DiagonalMatrix[RR];  Clear[\[Psi]]  usol = NDSolveValue[{I D[\[Psi][t], t] ==      Z.\[Psi][t], \[Psi][0] == \[Psi]ini}, \[Psi], {t, tini, tmax}]  Plot[usol[t], {t, tini, tmax}] 

Now, I´m trying to solve the same system but writing the equations

Clear["Global*"]  tini = 0;  tmax = 200;  usol = NDSolveValue[{I x1'[t] == x2[t], I x2'[t] == x1[t] + x3[t],      I x3'[t] == x2[t] + x4[t], I x4'[t] == x3[t] + x5[t],      I x5'[t] == x4[t], x1[0] == 0, x2[0] == 0, x3[0] == 0, x4[0] == 0,      x5[0] == 1}, {x1, x2, x3, x4, x5}, {t, tini, tmax}];  Plot[usol[t], {t, tini, tmax}] 

Why the second code doesn´t give me the same result if I write the same system?

## Ndsolve of random (periodic) initial conditions via BSplineFunction

Suppose we want to solve the heat equation for a random initial temperature field T(t,x,y). Therefore I define the initial conditions with the BSplineFunction which nicely incorporates periodic boundary conditions and try to solve it with NDSolve:

Clear["Global*"] Lx = 1; Ly = 1; eq = D[T[t, x, y], t] == Laplacian[T[t, x, y], {x, y}]; pbc1 = T[t, x, 0] == T[t, x, Ly]; pbc2 = T[t, 0, y] == T[t, Lx, y]; T0[x_, y_] =   BSplineFunction[RandomReal[1, {30, 30, 1}], SplineClosed -> True][x,    y] Plot3D[T0[x, y], {x, 0, Lx}, {y, 0, Ly}, PlotRange -> All,   AxesLabel -> {"x", "y", "T"}, ColorFunction -> "DarkRainbow",   Mesh -> All, MeshStyle -> Opacity[.2]] T0[0, 0] ic = T[0, x, y] == T0[x, y] ic2 = T[0, x, y] == 0.5*Exp[-4 ((x - Lx/2)^2 + (y - Ly/2)^2)]; Monitor[AbsoluteTiming[   sol = NDSolve[{eq, pbc1, pbc2, ic},      T, {t, 0, 1}, {x, 0, Lx}, {y, 0, Ly},      Method -> {"MethodOfLines",        "SpatialDiscretization" -> {"TensorProductGrid",          "MaxPoints" -> 25}},      EvaluationMonitor :> (currentTime =         Row[{"t = ", CForm[t]}])]], currentTime] 

The code unfortunately throws an error "Data … is not a rectangular tensor with dimensions …". I assume this results from the fact that T0 evaluated e.g. at grid point [0,0] returns a list {} and NDSolve probably needs just the float value.

I tried T0[0,0][[1]] which yields the expected value. But this does not work for the definition of the initial condition: ic = T[0, x, y] == T0[x, y][[1]], where x,y are not set yet. This returns just x, in reality it should first put in the values of x,y and only after that take list element 1 (if I am not mistaken).

Apart from that the code works fine which can be verified if you simply use ic2 instead of ic as initial condition in NDSolve.

Can someone help to get this running?

## Derivative of NDSolve output

I am trying to study the stability of a 2-dimensional discrete system (X, Y) by finding the Jacobian at the systems non-trivial equilibrium (X*, Y*). The functions that map the system from one iterate, n, to the next are given by the solutions to ordinary differential equations solved at time t = T

X(n+1) = f(x(T), y(T))

Y(n+1) = g(x(T), y(T))

where x(T) and y(T) depend on X(n) and Y(n).

Here is a MWE of the problem

sol[T_?NumericQ, a_, b_, c_, x0_?NumericQ, y0_?NumericQ] :=    NDSolve[{     x'[t] == x0 a - b *x[t]*y[t],     y'[t] == b*x[t]*y[t] - c y[t],     x[0] == 0, y[0] == y0},     {x, y}, {t, 0, T];  xmap[T_?NumericQ, a_, b_, c_, x0_?NumericQ, y0_?NumericQ] :=   x[T] /. sol[T, a, b, c, x0, y0] ymap[tt_?NumericQ, a_, b_, c_, h_, j_, x0_?NumericQ, y0_?NumericQ] :=   h (y[T] /. sol[T, a, b, c, x0, y0])/(1 +       j (y[T] /. sol[T, a, b, c, x0, y0])) 

where x0 and y0 are placeholders for X(n-1) and Y(n-1). Is it possible to take the partial derivatives of xmap and ymap w.r.t. x0 and y0? I tried this

D[xmap[5, 0.5, 0.25, 1, x0, 200], x0] /. x0 -> 100 

But it does not evaluate.

## NDSolve How do I insert another boundry condition

So I’m new to programming and I am trying to solve a very simple differential equation (i can do this by hand but the thing is i can’t program it), the equation is a'[t] == - 0.0118, with the boundary conditions a[0]==3 and I want to know a[180]=?, how do I insert this (a[180]=?) into the code, this is my short code.

Sol = NDSolve[{a'[t] == -0.0118 , a[0] == 3 }, a, {t, 0, 180}] Plot[Sol, {t, 0, 180}]  

It doesn’t show anything in the plot

Thanks a lot, any type of help will be appreciated

## Solving simultaneous ODEs with NDSolve causes unexpected singularity error (ndsz)

I’m trying to model the evaporation of a binary droplet in a square well by solving a couple of PDEs for its height profile and composition. I can manage it for a single liquid pretty well, but when I introduce the PDE for composition NDSolve quickly finds a singularity in the composition variable (importantly the singularity moves slightly if I change the number of points). I can’t seem to find any way to stabilise the system, can anyone help me? Is there something fundamental about NDSolve that I’m missing?

(I asked a similar question a while ago but I’ve restructured my code to make it easier to share here)

The equations are:

$$\frac{\partial h}{\partial t} = – C \frac{\partial}{\partial x}\Big[\frac{1}{3} \frac{\partial^3h}{\partial x^3}h^3 + \frac{B1}{2} \frac{\partial X}{\partial x}\Big] – (1-AX)$$ $$\frac{\partial X}{\partial t} = -\frac{C}{3} h^2 \frac{\partial^3h}{\partial x^3}\frac{\partial X}{\partial x} – \frac{hB1}{2} \Big(\frac{\partial X}{\partial x}\Big)^2 + \frac{A}{h}X(1-X)$$

with height $$h$$ and composition $$X$$.

And the boundary conditions are

$$h(t)=1$$ and $$X(t)=X(t=0)$$ at $$x=1$$, $$\frac{\partial h}{\partial x}=\frac{\partial X}{\partial x}=0$$ at $$x=0$$, and $$\frac{\partial^3h}{\partial x^3}=0$$ at $$x=0$$ and $$x=1$$.

I discretise all this using finite differences.

This is what I get for a single liquid droplet (when I suppress equation 2) and I expect something qualitatively similar for a binary droplet.

Here is a minimum replicable example:

    N1 = 100; dx = 1/N1; (*Discretise*)     C1  = 1; (*parameter C*)     a = 1.5; b = 1 - a; hInitial = Table[a + b i^2 dx^2, {i, 0, N1}]; (*Initial height profile*)     E1 = 1;(*evaporation on/off switch*)     Cc = C1/(24 * dx^3);     A1 = 0.01; B1 = 1;  (*Control constants in the equations above*)     CcX = (C1 B1)/(2 dx^2); CcX2 = C1/(6 dx^2);     X0 = 0.4; XInitial = Table[X0, {i, 0, N1}]; (*Initial composition ratio*)     Sc = 1; (*Composition on/off switch*)     TGap = 0.085;(*for spacing curves in the plot*)         dv[v_List] :=       With[{h = Take[v, Length[v]/2 - 1], hN = v[[(Length[v]/2)]],          X1 = v[[Length[v]/2 + 1 ;; Length[v] - 1]], XN = v[[-1]]},        With[{          dh1 =            ListCorrelate[{0, 0, 1 , 1, 0}, #] &, (*derivatives are discretised as finite differences.            ListCorrelate achieves this.*)           dh2 = ListCorrelate[{0, 1, 1, 0, 0}, #] &,          dh3 = ListCorrelate[{0, -1, 3, -3, 1}, #] &,          dh4 = ListCorrelate[{-1, 3, -3, 1, 0}, #] &,          hi = ListCorrelate[{0, 0, 1, 0, 0}, #] &,          dh5 = ListCorrelate[{-2, 1, 0, -1, 2}, #] &,            dx1 = ListCorrelate[{0, -1, 1}, #] &,          dx2 = ListCorrelate[{-1, 1, 0}, #] &,          Xi = ListCorrelate[{0, 1, 0}, #] &          },         Flatten@{             -Cc*(dh1[#1]^3*dh3[#1] - dh2[#1]^3*dh4[#1]) -               CcX*(dh1[#1]^2*dx1[#2] - dh2[#1]^2*dx2[#2]) -              E1*(1 - A1 Xi[#2]), (*ODE for heights*)              0,(*height is pinned at edge, (1, 1)*)              Sc*(-CcX2* hi[#1]^2*dh5[#1]*dx2[#2] -              CcX/C1* hi[#1]*dx2[#2]^2 +              A1 *(hi[#1])^-1*Xi[#2]*(1 - Xi[#2])), (*ODE for composition*)                      0 (*composition ratio is fixed at edge, XN = 0.4*)             } &[          Flatten@Join[{h[[3]]}, {h[[2]]}, {h}, {1}, {3 - 3 h[[-1]] + h[[-2]]}],                       Flatten@Join[{X1[[2]]}, {X1}, {XN}]]]                (*this contains the boundary conditions*)         ];      v0 = Join[hInitial, XInitial]; (*Initial conditions*)     system2d =       NDSolve[{v'[t] == dv[v[t]], v[0] ==  v0,          WhenEvent[          Min@Table[             Take[values[[tt]], Length[values[[tt]]]/2],               {tt, 1,  Length[values]}] < 0,           "StopIntegration"](*WhenEvent to stop integration at touchdown*)},          v, {t, 0, 2}][[1, 1, 2]];       Needs@"DifferentialEquationsInterpolatingFunctionAnatomy";     values = InterpolatingFunctionValuesOnGrid@system2d;      times = Flatten@InterpolatingFunctionGrid@system2d;      T1 = 0; timesteps = {}; For[i = 1, i < Length[times], i++,       If[times[[i]] > T1, {T1 = T1 + TGap,         timesteps = Append[timesteps, {i, times[[i]]}]}]];     finaltime = {{-1}}; firstandfinaltime = {{1}, {-1}};      heights =        Table[Take[values[[tt]], Length[values[[tt]]]/2],          {tt, 1, Length[values]}];      comps = Table[        Take[values[[tt]], -(Length[values[[tt]]]/2)],           {tt, 1, Length[values]}];      Show[Table[       ListPlot[Table[{i*dx, heights[[ttt[[1]], i + 1]]}, {i, 0, N1}],         PlotRange -> {{0, 1}, {0, 1.05*a}}, Joined -> True,         AxesLabel -> {x, h}], {ttt, finaltime}]]      Show[Table[       ListPlot[Table[{i*dx, comps[[ttt[[1]], i + 1]]}, {i, 0, N1}],         PlotRange -> {{0, 1}, {0, 1}}, Joined -> True,        AxesLabel ->{x, X}], {ttt, finaltime}]] 

## Why both DSolve and NDSolve are unable to solve a second-order differential equation?

I am trying to solve a recurrence relation using generating functions method. After some long calculations, I have arrived to this second-order differential equation: $$0.5 x^5 y”(x)+(2x^4+x^3)y'(x)+\left(x^3+x^2+x-1\right)y(x)+1=0$$

and these conditions: $$y(0)=1, y'(0)=1$$. $$y(x)$$ is the function that needs to be expanded as Taylor Series at $$x=0$$ to obtain the sequence from the coefficients. However, when I try to solve it both using DSolve and NDSolve, I have no luck. With DSolve it just returns the request itself:

$$\text{DSolve}\left[\left\{0.5 x^5 y”(x)+(2. x+1) x^3 y'(x)+\left(1. x^3+x^2+x-1\right)y(x)+1=0,y(0)=1,y'(0)=1\right\},y,x\right]$$

And with NDSolve I just receive errors and no equation:

Power::infy: Infinite expression 1/0.^5 encountered. Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.. 

$$\text{NDSolve}\left[\left\{0.5 x^5 y”(x)+(2. x+1) x^3 y'(x)+\left(1. x^3+x^2+x-1\right)y(x)+1=0,y(0)=1,y'(0)=1\right\},y,\{x,0,1\}\right]$$

How could I resolve this problem?

## When does NDSolve issue zero step size error NDSolve::ndsz?

I want to detect directly when step size becomes "effectively" zero. The below example from the documentation throws the error message as expected:

s = {}; NDSolve[{(2 - f[x]) f'[x] == f[x], f[0] == 1}, f, {x, 0, 5}, StepMonitor :> AppendTo[s, x]];  NDSolve::ndsz: At x == 0.3862940268757776, step size is effectively zero; singularity or stiff system suspected. 

The below code indicates that none of the actual steps taken has zero length.

AnyTrue[Differences@s, PossibleZeroQ]  (* False *) 

How does NDSolve decide that the step size is zero? I can of course capture the NDSolveValue::ndsz error, but I want to know when exactly (depending on what parameters) the error is issued. In some extreme cases, NDSolve may generate an InterpolatingFunction solution that has practically zero domain length (but not according to PossibleZeroQ).

## Improving mesh and NDSolve solution convergence

I have developed the code below to solve two PDEs; first mu[x,y] is solved for, then the results of mu are used to solve for phi[x,y]. The code works and converges on a solution as is, however, I would like to decrease the size of a, b, and d even further. To accurately represent the physical process I am trying to simulate, a, b, and d would need to be ~100-1000x smaller. If I make them smaller, I don’t believe the solution has actually converged because the values for phi along the right boundary change significantly with a change in mesh size (i.e. if I make them smaller and the code below produces a value of phi=-0.764 at the midpoint between y2 and y3 along the right boundary, a change in size1 to 10^-17 and size2 to 10^-15, changes that value of phi to -0.763, and a change in size2 to 10^-16 changes that value again to -0.860), but I cannot make the mesh size any smaller without Mathematica crashing.

Are there any better ways to create the mesh that would be less computationally taxing and allow it to be more refined in the regions of interest? Or are there any ways to make the code in general less computationally expensive so that I can further refine the mesh?

ClearAll["Global*"] Needs["NDSolveFEM"] (* 1) Define Constants*) e = 1.60217662*10^-19; F = 96485; kb = 1.381*10^-23; sigi = 18; sigini = 0; sigeni = 2*10^6; T = 1000; n = -0.02; c = 1;  pH2 = 0.2; pH2O = 1 - pH2; pO2 = 1.52*^-19; l = 10*10^-6; a = 100*10^-7; b = 50*10^-7; d = 300*10^-7; y1 = 0.01; y2 = 0.5*y1; y3 = y2 + a; y4 = y3 + d; y5 = y4 + b; mu1 = 0; mu2 = -5.98392*^-19; phi1 = 0;  (* 2) Create mesh*) m = 0.1*l; size1 = 10^-16; size2 = 10^-15; size3 = 10^-7; mrf = With[{rmf =       RegionMember[       Region@RegionUnion[Disk[{l, y2}, m], Disk[{l, y3}, m],          Disk[{l, y4}, m], Disk[{l, y5}, m]]]},     Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices];      Which[rmf[{x, y}],        area > size1, (0 <= x <= l && y2 - l <= y <= y2 + l),        area > size2, (0 <= x <= l && y3 - l <= y <= y3 + l),        area > size2, (0 <= x <= l && y4 - l <= y <= y4 + l),        area > size2, (0 <= x <= l && y5 - l <= y <= y5 + l),        area > size2, True, area > size3]]]]; mesh = DiscretizeRegion[Rectangle[{0, 0}, {l, y1}],     MeshRefinementFunction -> mrf];  (* 3) Solve for mu*) bcmu = {DirichletCondition[mu[x, y] == mu1, (x == 0 && 0 < y < y1)],    DirichletCondition[     mu[x, y] ==       mu2, (x == l && y2 <=  y <=  y3) || (x == l && y4 <= y <= y5)]}; solmu = NDSolve[{Laplacian[mu[x, y], {x, y}] ==       0 + NeumannValue[0, y == 0 || y == y1 ||         (x == l && 0 <= y < y2) || (x == l &&            y3 < y < y4) || (x == l && y5 < y < y1)], bcmu},     mu, {x, y} \[Element] mesh, WorkingPrecision -> 50];  (* 4) Solve for electronic conductivity everywhere*) pO2data = Exp[(mu[x, y] /. solmu)/kb/T]; sige0 = 2.77*10^-7; sigedata = Piecewise[{{sige0*pO2data^(-1/4), 0 <= x <= l - m},     {sige0*pO2data^(-1/4), (l - m < x <= l && 0 <= y < y2)},     {(sigeni - sige0*(pO2data /. x -> l - m)^(-1/4))/m*(x - (l - m)) +        sige0*(pO2data /. x -> l - m)^(-1/4), (l - m < x <= l &&         y2 <=  y <= y3)},     {sige0*pO2data^(-1/4), (l - m < x <= l && y3 < y < y4)},     {(sigeni - sige0*(pO2data /. x -> l - m)^(-1/4))/m*(x - (l - m)) +        sige0*(pO2data /. x -> l - m)^(-1/4), (l - m < x <= l &&         y4 <= y <= y5)},     {sige0*pO2data^(-1/4), (l - m < x <= l && y5 < y <= y1)}}];  (* 5) Solve for phi*) Irxn = -(2*F)*(c*pO2^n ); A = (Irxn - sigi/(4*e)*(D[mu[x, y] /. solmu, x] /. x -> l))/(-sigi); B = sigi/(4*e)*(D[mu[x, y] /. solmu, x] /.        x -> l)/(sigi + sigedata /. x -> l - m); bcphi = DirichletCondition[phi[x, y] == phi1, (x == 0 && 0 < y < y1)]; solphi = NDSolve[{Laplacian[phi[x, y], {x, y}] ==       0 + NeumannValue[0,         y == 0 ||          y == y1 || (x == l && 0 <= y < y2) || (x == l &&            y3 < y < y4) || (x == l && y5 < y < y1)] +        NeumannValue[-A[[1]], (x == l && y2 <= y <= y3)] +        NeumannValue[-B[[1]], (x == l && y4 <= y <= y5)], bcphi},     phi, {x, y} \[Element] mesh, WorkingPrecision -> 50];  (* 6) Print values to check for convergence*) P[x_, y_] := phi[x, y] /. solphi; P[l, (y3 - y2)/2 + y2] P[l, (y5 - y4)/2 + y4] 

## Exiting CrossSlidingDiscontinuity in NDSolve to follow equilibrium curve

I am trying to figure out how to get the solution curve to an NDSolve to slide along once it reaches a boundary and then to exit the CrossSlidingDiscontinuity once it reaches an equilbrium curve and follow the equilibrium curve back within a set boundary.

If the x[t] solution curve reaches the boundary edge on this interval [0,1] I needed it to slide along this boundary, hence I used CrossSlidingDiscontinuity. However once, the solution reaches the yellow equilibrium sine curve, I need it to exit the sliding discontinuity and follow the equilibrium curve. I tried using EventLocator (in the code below). It manages to fall back within [0,1] but does not follow the equilibrium sine curve. I attached the image of the plot below.

I also tried using two WhenEvents (one for when x[t] reaches 1 triggering the sliding and one where x[t] reaches the sine curve triggering to follow the sine curve) but only the sliding WhenEvent was triggered.

Here is my code:

normalDE2[x_, t_] := -(x + 1/2)*(x - 1/2)*(x - (1 + .1*Sin[t]));  testDE3[x_?NumberQ, t_] :=  If[0 < x < 1, normalDE2[x, t], -1*normalDE2[x, t]];  sol4 = NDSolve[{x'[t] == testDE3[x[t], t], x[0] == .75,  WhenEvent[{x[t] == 1, x[t] == 0},  {Print[t],"CrossSlidingDiscontinuity"}]},  x, {t, 0, 10}, Method -> {"EventLocator", "Event" -> {x[t] - (1 +.1*Sin[t])}, "EventAction" :> {{x[t] -> (1 + .1*Sin[t])}}}];  Plot[{x[t] /. sol4, (1 + .1*Sin[t])}, {t, 0, 10},  Frame -> True, GridLines -> Automatic]  

Plot Image