Solving PDE with prescribed rectangular mesh

Question: I want to solve a PDE using NDEigensystem. The region is a rectangle Rectangle[{0,-70}, {20, 70}]. How can I make a mesh that satisfies the following two conditions?

(1) The mesh is finer near $ y \approx 0$ .

(2) The line $ y=0$ should be the boundary of the mesh.

Motivation: The PDE that I want to solve can be written in the following command: First, I set the mesh as below:

mesh = ToElementMesh[Rectangle[{0, -70}, {20, 70}],     MeshRefinementFunction ->     Function[{vertices, area},      Block[ {x, y}, {x, y} = Mean[vertices];       If[-10 < y < 10, area > 0.1, area > 10]]] ] 

and then write down the PDE:

NDEigensystem[{{3*psi3[x, y] - 10*psi2[x, y]*Sign[y] - Derivative[0, 1][psi2][x, y] -      I*Derivative[1, 0][psi2][x, y], 3*psi4[x, y] - 10*psi1[x, y]*Sign[y] + Derivative[0, 1][psi1][x, y] -      I*Derivative[1, 0][psi1][x, y], 3*psi1[x, y] - 10*psi4[x, y]*Sign[y] - Derivative[0, 1][psi4][x, y] -      I*Derivative[1, 0][psi4][x, y], 3*psi2[x, y] - 10*psi3[x, y]*Sign[y] + Derivative[0, 1][psi3][x, y] -      I*Derivative[1, 0][psi3][x, y]}, PeriodicBoundaryCondition[psi1[x, y], x == 0,     TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], PeriodicBoundaryCondition[psi2[x, y], x == 0,     TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], PeriodicBoundaryCondition[psi3[x, y], x == 0,     TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], PeriodicBoundaryCondition[psi4[x, y], x == 0,     TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], DirichletCondition[psi1[x, y] == 0,     0 < x < 20 && (y == 70 || y == -70)], DirichletCondition[psi2[x, y] == 0,     0 < x < 20 && (y == 70 || y == -70)], DirichletCondition[psi3[x, y] == 0,     0 < x < 20 && (y == 70 || y == -70)], DirichletCondition[psi4[x, y] == 0,     0 < x < 20 && (y == 70 || y == -70)]}, {psi1, psi2, psi3, psi4}, Element[{x, y}, mesh], 10] 

The single important fact about the PDE is that the differential operator contains $ \mathrm{sgn}(y)$ , hence it is not continuous at $ y=0$ . This is why I want the mesh boundary to include $ y=0$ . The mesh I created is visualized as follows:

enter image description here

As wanted, the mesh is finer. However, $ y=0$ is not the boundary of the mesh.

The result of one solution obtained by solving the PDE is plotted as below:

enter image description here

In the above plot, the horizontal direction is the $ x$ -direction. Note that the solution is localized near $ y\approx 0$ . This is why I want to make the mesh finer near $ y\approx 0$ . I don’t believe the result is the good representation of the solution, since the graph looks jagged, which I don’t expect since the differential operator has a translation symmetry in $ x$ -direction.

Solving a complicated time-dependent function

I want to solve the following time-dependent system of equations:

Clear["Global`*"]; q = (((1602176634)/10^9))*10^(-19); k = (((1380649)/10^6))*10^(-23); \[Eta] = 3/2; Td = 20 + (5463/20); R1 = 10; Vi = Piecewise[{{u*Sin[\[Omega]*t],       0 <= t <= ((2*Pi)/\[Omega])/2}, {0, t > ((2*Pi)/\[Omega])/2}}]; u = 230*Sqrt[2]; \[Omega] = 2*Pi*50; Is = 5*10^(-9); FullSimplify[  Solve[{I1 == Is*(Exp[(q*(Vi - V1))/(\[Eta]*k*Td)] - 1),     I1 == Is*(Exp[(q*(V1 - V2))/(\[Eta]*k*Td)] - 1), I1 == V2/R1}, {I1,     V1, V2}]] 

But it spits out nothing, just the simplified version of the input. Is there a way to solve for the unknowns $ I_1$ , $ V_1$ and $ V_2$ ?

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. Droplet curves when the composition equation is suppressed

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@"DifferentialEquations`InterpolatingFunctionAnatomy`";     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}]] 

Thank you for any help you can offer me!

Numerically solving 3D Maxwell equations with NDEigensystem

I am trying to get electric $ \vec{E}$ and magnetic $ \vec{B}$ fields in a cylindrical cavity with a dielectric, as in the following Figure. Cavity (pink) with a dielectric inside (blue).

Both the cavity (pink) and dielectric (blue, with dielectric constant \epsilon_r) are cylindrical and share the main axis. The cavity is assumed to be conducting, such as the field at its surface has to be null.

By resorting to the vector potential $ \vec{A}$ , and using the generalized Coulomb Gauge transformation such that:

$ \vec{E}(t,\vec{r})=-\partial_t \vec{A}(t,\vec{r}),$

$ \vec{B}(t,\vec{r})=\vec{\nabla} \times \vec{A}(t,\vec{r}),$

$ \vec{\nabla} \cdot \left[ \epsilon_r(\vec{r}) \vec{A}(t,\vec{r}) \right]= 0,$

I got the system of differential equations that $ \vec{A}$ has to satisfy: \begin{equation} \vec{\nabla}^2 \vec{A}(t,\vec{r}) – \frac{1 + \epsilon_r(\vec{r})}{c^2}\partial_t^2 \vec{A}(t,\vec{r})=0. \end{equation} Here, $ \vec{r} = (x,y,z)$ .

I would like to have the eigenfrequencies and (spatial) eigenfunctions of this last operator. As far as I have understood, NDEigensystem with DirichletCondition is what I have to use; with the first output being the frequencies of the modes and the second output the spatial envelopes of the eigenfunctions.

I tried that, and failed (quite miserably). Not only the eigenfrequencies are imaginary and the spatial envelopes of the eigenfunctions do have an imaginary part$ ^{*}$ . If I plug the outcome of NDEigensystem into the differential equations, I find that these equations are not even satisfied. I am sure that it is me failing somewhere, but after long time spent trying, I am starting being really frustrated. My code is the following:

{vals, funs} = NDEigensystem[{EqSt[t, x, y, z, z1, z2, e, c] == {0, 0, 0}, BndCnd}, {Ax[t, x, y, z], Ay[t, x, y, z], Az[t, x, y, z]}, t, {x, y, z} \[Element] Cylinder[{{0, 0, 0}, {0, 0, d}}, r], 16, Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}}]; 

where EqSt[t, x, y, z, z1, z2, e, c] is the system of differential equations

EqSt[t_, x_, y_, z_, z1_, z2_, er_, c_] := Laplacian[{Ax[t, x, y, z], Ay[t, x, y, z], Az[t, x, y, z]}, {x, y, z}] - (1 + fer[z, z1, z2, er])/c^2 {D[Ax[t, x, y, z], {t, 2}], D[Ay[t, x, y, z], {t, 2}], D[Az[t, x, y, z], {t, 2}]}; 

and BndCnd the boundary conditions

BndCnd = {DirichletCondition[Ax[t, x, y, z] == 0, True], DirichletCondition[Ay[t, x, y, z] == ,True], DirichletCondition[Az[t, x, y, z] == 0, True]}; 

Finally, the dielectric function fer[z, z1, z2, er] that I use to mimic the dielectric is the following:

fer[z_, z1_, z2_, e_] := e (HeavisideTheta[z - z1] - HeavisideTheta[z - z2]); 

which is a step function in the z-th coordinate (= the axis of both cylinders) having value e betweem z1 and z2 and zero elsewhere.

I have tried different methods and different values of "MaxCellMeasure" (there is little improvement, but the error which I obtain by plugging the solution into the differential equations is above one!!!). Do you have any idea what is wrong here?


I have few "bonus questions" here.

-First, I tried increasing the mesh resolution for getting better outcomes, but all the times the eigenfunctions look quite bad, not smooth at all. I guess that a 3D mesh is quite demanding, but I know that the most difficult part is going to be at the dielectric. How can I tell NDEigensystem to have a finer mesh there?

-Second, there is a way to change the normalization of the eigenfunctions? I have read that NDEigensystem spits out eigenfunctions $ \vec{\phi}_i$ s.t.: $ \int \vec{\phi}_i^* \cdot \vec{\phi}_j d \vec{r} = \delta_i,j$ . I would change that to $ \int \epsilon_r(\vec{r}) \vec{\phi}_i^* \cdot \vec{\phi}_j d \vec{r} = \delta_i,j$ , would this be possible?

-Finally, I guess that it might be easier for NDEigensolver to remove the time dependence from the differential equations. This can be easily done by assuming that $ \vec{A}(t,\vec{r}) \rightarrow \vec{A}(\vec{r})e^{-i \omega t})$ . In this case, one can rewrite the operator equation above as: \begin{equation} \vec{\nabla}^2 \vec{A}(\vec{r}) + \frac{\omega^{2}}{c^2}[1 + \epsilon_r(\vec{r})] \vec{A}(\vec{r})=0. \end{equation} However, I do not know how to tell NDEigensystem that $ \omega$ is not a parameter, but should be found with the system constraints… There is a way to remove the time dependence from the equation to be given in NDEigensystem?


$ ^*$ I guess this is a stupid problem; if I change the sign of the Laplacian I get real eigenfrequencies and eigenfunctions… But I am quite confident that the sign in the above equation is correct. There is something that I am missing here?

Solving trigonometric equations with two variables in fixed range?

I am trying to use ‘Solve’ function to solve two trigonometric equations with two variables a1 and a2 (with range 0,2Pi). I wrote the code in the following:

  Vets={Cos[2a1](I+Cos[2a2])+Sin[2a1]Sin[2a2], (I-Cos[2a2])Sin[2a1]+Cos[2a1]Sin[2a2]};   Solve[Vets[[1]] == 1 && Vets[[2]] == 0 && 0<=a1<=2Pi && 0<=a2<=2Pi, {a1,a2}]  

but it gives errors:

Solve::nsmet: This system cannot be solved with the methods available to Solve.

I looked up the documents for Solve function and it should be no problem with the above code. But it doesn’t do anything and give errors, I don’t understand why?

It will be super cool if anyone could answer this question, thank you very much in advance!

Solving Kronecker-structured linear equations

I need to approximately solve the following underdetermined system of $ n$ linear equations

$ $ y_i=a_i^T X b_i$ $

Where $ X$ is $ d\times d$ unknown matrix, $ a_i$ and $ b_i$ are given vectors and $ n=d$ . Is there a way to do this much faster than vectorizing both sides of each equation and calling LinearSolve?

LinearSolve approach destroys the structure, it reduces to solving a system with $ d^3$ coefficients instead of original one with $ 2d^2$ coefficients.

Below is an example of this approach, it’s too slow for my application where $ d=1000$ . On that scale, backward is 2000x slower than forward, with majority of the time spent in LinearSolve. I was hoping for 100x slower since that seems like typical LinearSolve overhead for unstructured systems on that scale.

n = d = 50; {a, b} = Table[RandomReal[{-1, 1}, {n, d}], {2}]; X = RandomReal[{-1, 1}, {d, d}]; forward[X_] := MapThread[#1.X.#2 &, {a, b}]; y = forward[X];  backward[Y_] :=    With[{mat = MapThread[Flatten@Outer[Times, #1, #2] &, {a, b}]},    x = LinearSolve[mat, y];    ArrayReshape[x, {d, d}]    ]; Print["forward time is ", forward[X]; // Timing // First]; {timing, error} = Timing[Norm[forward[backward[y]] - y, Infinity]]; Print["backward time is ", timing]; Print["error is ", error] 

Solving the heat equation using Laplace Transforms

I am trying to solve the 1-D heat equation using Laplace Transform theory. The equation is as follows. I don’t have the capability to write the symbols so I will write it out.

                     partial u/partial t = 2(partial squared u/ partial x squared) -x    boundary conditions are partial u/partial x(0,t)=1, partial u/partial x(2,t)=beta. 

The problem asks the following: (a). For what value of beta does there exist a steady-state solution? (b). if the initial temperature is uniform such that u(x,0)=5 and beta takes the value suggested by the answer to part (a), derive the equilibrium temperature distribution.

I was able to get an equation that looks like U(x,s)=c e^(s/2)^1/2 -(1/s)((x/s)-u(x,0)). But I am not sure how to go from here to solve for beta using the boundary conditions. I need some assistance from someone.

Solving a system of differential equations whose one of the coefficients is imported data

Suppose we have a coupled system of differential equations: \begin{equation} \frac{db}{dt}=(- \gamma_b -i\omega_b)b-i\frac{g}{2}p;\quad \frac{dp}{dt}=i\frac{g}{2}\Delta N(t) b-(\gamma_a+\gamma_b+2iJ)p. \end{equation} If $ \Delta N$ was fixed, the solution of the system would be like \begin{equation} \begin{pmatrix} b(t)\ p(t) \end{pmatrix}=\begin{pmatrix} a_{11}&a_{12}\ a_{21}&a_{22} \end{pmatrix}\begin{pmatrix} b(0)\ p(0) \end{pmatrix} \end{equation} Using the following code, I have found a $ 2\times 2$ matrix (called sol) whose entries are $ a_{ij}$ in the above equation:

rb=630;wb=75*10^6;g=0.63;ra=2.6*10^6;rm=3.6*10^6;J=6.3*10^7;DeltaN=0.164*10^5; m ={{-rb-I wb,-I g/2},{I g DeltaN/2,-(ra+rm+2 I J)}}; eigvec = Eigenvectors[m] // Transpose // Simplify; eigval = Eigenvalues[m] // Simplify; inv = Inverse[eigvec] // Simplify; v1 = eigval[[1]]; v2 = eigval[[2]]; sol = eigvec.{{E^(v1 t), 0}, {0, E^(v2 t) }}.inv; 

If we suppose that $ p(0)=0$ , then one can easily plot $ |b(t)/b(0)|^2$ : simply plot $ a_{11}(t)$ . But the problem is that $ \Delta N$ is not fixed. It is a $ N\times 1$ matrix which I have obtained from another code written with Fortran and its type is data.txt. The elements of this file are calculated by assuming the time interval between each one is $ 0.001$ . That is, for $ t=0.001$ we have $ \Delta N_1$ , for $ t=0.002$ we have $ \Delta N_2$ , etc. But the time intervals are not included in the txt file.

One way that comes to my mind is this: Assuming we know the analytical form of solfor a fixed $ \Delta N$ , we set time, i.g., equal to $ 0.001$ and then substitute the first row of the txt file (I call it $ \Delta N_1$ ) into sol and find $ a_{11}$ . Then we raise time to $ 0.002$ , substitute $ \Delta N_1$ into sol, find $ a_{11}$ , and repeat the procedure to the last row of the txt file.

Now the question is this: how can I import the txt file to the code and do the procedure that I explained above to get some data like $ \{\{0.001,a11(0.001)\},\{0.002,a11(0.002)\},….\}$ where the first elements are time intervals and the second ones are $ a_{ij}$ corresponding to that particular time?

I had asked a similar question here enter link description here, but in that problem I did not have an external file with txt format.

I could not upload my txt file, so I write the first 10 elements if necessary:

0.164E+05

0.655E+05

0.146E+06

0.258E+06

0.400E+06

0.572E+06

0.776E+06

0.101E+07

0.129E+07

0.159E+07

In what cases is solving Binary Linear Program easy (i.e. **P** complexity)? I’m looking at scheduling problems in particular

In what cases is solving Binary Linear Program easy (i.e. P complexity)?

The reason I’m asking is to understand if I can reformulate a scheduling problem I’m currently working on in such a way to guarantee finding the global optimum within reasonable time, so any advice in that direction is most welcome.

I was under the impression that when solving a scheduling problem, where a variable value of 1 represents that a particular (timeslot x person) pair is part of the schedule, if the result contains non-integers, that means that there exist multiple valid schedules, and the result is a linear combination of such schedules; to obtain a valid integer solution, one simply needs to re-run the algorithm from the current solution, with an additional constraint for one of the real-valued variables equal to either 0 or 1.

Am I mistaken in this understanding? Is there a particular subset of (scheduling) problems where this would be a valid strategy? Any papers / textbook chapter suggestions are most welcome also.