14

Update (Steady-State Solution) I think the fundamental issue is that you are over constraining your system. Whether you are solving the "heat equation" or not, your operator has the same form of the heat equation as shown below: $$\rho {{\hat C}_p}\frac{{\partial T}}{{\partial t}} + \nabla \cdot {\mathbf{q}} = 0$$ If the flux, $\mathbf{q}$, needs ...


12

DSolve can handle this. Clear[y]; y[x_, e_] = y[x] /. DSolve[{ e y''[x] + y'[x] x == -e Pi^2 Cos[Pi x] - Pi x Sin[Pi x], y[-1] == -2, y[1] == 0}, y[x], x][[1]] // Simplify Manipulate[ Plot[y[x, e], {x, -1, 1}], {{e, 0.01}, 0.0001, 0.1, Appearance -> "Labeled"}]


12

I think here is another way to do it. For this we use the low level FEM functions. First we have a utility function that converts a PDE specification to a discrete version. Needs["NDSolve`FEM`"] PDEtoMatrix[{pde_, \[CapitalGamma]___}, u_, r__, opts : OptionsPattern[NDSolve`ProcessEquations]] := Module[{ndstate, feData, sd, vd, bcData, methodData, ...


12

The "Polar" coordinates programmed into Mathematica implicitly assume that the $\theta$ coordinate runs between -π and π, not between 0 and 2π as your initial ParametricPlot assumed: CoordinateChartData["Polar", "CoordinateRangeAssumptions"] (* #1[[1]] > 0 && -\[Pi] < #1[[2]] <= \[Pi] & *) In particular, ...


11

As @Henrik Schumacher points out, you have a very high aspect ratio (1000:1) domain. It is always conducive to conduct a dimensional analysis of your system. In the OP case, the dimensional analysis would show that the problem is essentially 1D in the $x$ direction. I will use the subscript $d$ to indicate that the variable has dimensions. First, we can ...


11

Two issues here. First of all, you've chosen 100 to approximate Infinity, which is way too large in this case. Something like 5 is OK: p = 0.2; inf= 5; boundaries = {-r + 1/2 (Sqrt[2] Sqrt[Cos[2 θ] (1 - p)^2 + 2 p + 1 - p^2] + 2 Cos[θ] (1 - p)), r - inf, -θ + Pi/2, θ - Pi}; << NDSolve`FEM` Ω = ToElementMesh@ImplicitRegion[And @@ (# <=...


10

As noted by @bbgodfrey, the "shooting" algorithms that Mathematica tends to use are not well-adapted to this particular equation. Better would be some kind of relaxation method, which is what Mathematica uses (I think) for solving PDEs on a mesh. And an ODE is just a PDE in one dimension, so let's try solving this equation on a one-dimensional mesh: Needs[...


10

Periodic potentials and Bloch waves This is a completely different approach, using functionality that was already present in Mathematica version 8. I'm posting this because it's quite robust and also allows me to demonstrate how Bloch's theorem can be used here to get a more general class of non-periodic solutions in a periodic potential. This is equivalent ...


10

Updated Background This is a more detailed response to a comment from the OP about the Neumann Value specification. The following is how I like to think of it. I was initially confused on how to setup Neumann values. After I read the Element Mesh Generation Tutorial, things got simpler for me because I could make the Mathematica workflow similar to other ...


9

DSolve is unable to solve this problem without assistance. Begin by solving the PDE only, s = DSolveValue[D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] == 4, u[x, y], {x, y}] /. {C[1] -> c1, C[2] -> c2} (* 2 x^2 + c1[I x + y] + c2[-I x + y] *) and then apply the two auxiliary conditions. b1 = Simplify[(s /. y -> x) == 2 x^2] (* c1[(1 + I) x] + c2[...


9

Solution in the case when the Neumann condition is given at all boundaries where possible shellID = 0.6; tubeLength = 6; rodOD = 0.03; rodCenters = shellID {{+0.25, -0.25}, {+0.25, +0.25}, {-0.25, -0.25}, {-0.25, \ +0.25}}; rodSection = RegionUnion @@ (Disk[#, rodOD] & /@ rodCenters); tubeSection = RegionDifference[Disk[{0, 0}, shellID], rodSection]; ...


9

Using DSolve V 12.1 can solve this exactly. ClearAll[u, x, y]; pde = Laplacian[u[x, y], {x, y}] == 0; L0 = 100; a = 1/10; b = 6/10; leftSide = Derivative[1, 0][u][0, y] == 0; rightSide = u[L0, y] == 0; bottomSide = Derivative[0, 1][u][x, 0] == a*u[x, 0]; topSide = u[x, L0] == b; bc = {leftSide, rightSide, bottomSide, topSide}; sol = DSolve[{pde, bc}, u[x,...


9

Linear PDEs typically are solved by the method of characteristics. For the PDE in the question, the ODEs of the characteristics are {x2'[s] == Sin[x1[s]] + x2[s], x1'[s] == x2[s], h'[s] == h[s]} Attempting to solve these ODEs with DSolve is unsuccessful, returning the Solve::ifun error message. That this is the same error cited in the question is not ...


8

To see why NDSolve has difficulty with this problem for very small e, consider that NDSolve solves this two-point boundary value problem by some form of shooting. In other words, it varies y'[min] until one is found that yields the desired y[max]. However, as e becomes very small, the sensitivity of y[max] to y'[min] becomes great, because the differential ...


8

In the apparent absence of a direct solution using NDEigensystem to the periodic problem posed in the question, ParametricNDSolveValue can be used. V[θ_] := Cos[θ]; s = ParametricNDSolveValue[{ϕ''[θ] + (V''[θ] - V'[θ]^2 - a) ϕ[θ] == 0, ϕ'[Pi/2] == 1, ϕ[0] == ϕ[2 Pi]}, ϕ, {θ, 0, 2 Pi}, {a}]; (The condition ϕ'[Pi/2] == 1 can, in principle, be applied ...


8

This does not get stuck, it returns right away ClearAll[x,y] z=1/1000; DSolve[{1+x-x^2==D[2/(y[x]*(((1+ z x)/(2*y[x]))^2+1)),{x}],y[0]==0},y[x],x] But your input makes it stuck ClearAll[x,y] z=0.001; DSolve[{1+x-x^2==D[2/(y[x]*(((1+ z x)/(2*y[x]))^2+1)),{x}],y[0]==0},y[x],x] As a general rules, use exact input with symbolic functions like DSolve, always ...


8

You can solve your problem by introducing a second ode (defining the antiderivative): sol = DSolve[{y''[x] == a y[x], Y'[x] == y[x], y'[L] == 0, Y[L] == 1}, {y, Y}, x][[1]] {y -> Function[{x}, (E^(-Sqrt[a] x) (E^(2 Sqrt[a] L) + E^(2 Sqrt[a] x)) C[1])/(1 + E^(2 Sqrt[a] L))], Y -> Function[{x}, (E^(-Sqrt[a] x) (Sqrt[a] E^(...


8

There is a simple sign error in the set up I would think; the left NeumannValues needs a negative sign: sol = NDSolveValue[{D[u[t, x], {t, 1}] - D[u[t, x], {x, 2}] == NeumannValue[-u[t, x], x == 0] + NeumannValue[1, x == 10], u[0, x] == 1 + x}, u, {x, 0, 10}, {t, 0, 10}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"...


8

If I understand the question right, you could use boundary markers like so: Needs["NDSolve`FEM`"] rego = RegionDifference[Disk[{0, 0}, 10], Disk[{0, 0}, 4]]; mesh = ToElementMesh[rego]; (*mesh["Wireframe"]*) The ToElementMesh auto generated markers. You can inspect them: groups = mesh["BoundaryElementMarkerUnion"] {1, 2} Visualize the boundary elements ...


8

This question is particularly interesting to me, and I have a package that may be helpful to you here. This particular equation is: Fourth order Linear Inhomogeneous in the independent variable Contains Robin-like boundary conditions Includes the eigenvalue in the boundary conditions I don't think that NDEigensystem can handle fourth-order systems, or ...


8

Too long for a comment. An easy way to generate a high quality mesh is to replace the Implicitegion with Cubuid and make use of the OpenCascade boundary mesh generator: Needs["NDSolve`FEM`"] (*a=ImplicitRegion[True,{{x,-1,1},{y,-1,1},{z,0,1}}];*) a = Cuboid[{-1, -1, 0}, {1, 1, 1}]; b = Cylinder[{{0, 0, -1/5}, {0, 0, 0}}, (650/1000)/2]; c = Cylinder[{{1, 1, ...


8

In a previous answer 240190, I showed how one could use anisotropic meshing to add a DirichletCondition at "infinity" for a 1D problem. In this answer, I shall extend the technique to a 2D problem. Geometry description In many FEM software packages, problems with spherical symmetry can be posed as an axisymmetric problem. Since it is easier for me ...


7

Here is a partial solution; especially to the issues of @bbgodfrey 's answer. Let's look at the fourth-order ODE: s1 = Flatten@ NDSolveValue[{ps''''[x] == 0, ps[0] == 0, ps'[0] == 0, ps[1] == 1, ps'[1] == 1}, {ps[x]}, {x, 0, 1}]; Plot[s1, {x, 0, 1}] This fourth-order ODE can be re-written as a system of two coupled second-order ODEs. s2 = Flatten@...


7

This problem can be solved symbolically, although perhaps not with DEigensystem. Instead, begin with DSolve. s = DSolveValue[{2 x y''''[x] + 4 y'''[x] == lamda y''[x], y[0] == 0, y'[1/2] == 0, y''[1/2] == lamda y[1/2]}, y[x], x] // FullSimplify (* (16 C[2] (BesselI[0, Sqrt[lamda]] (-1 + 2 x + Sqrt[2] Sqrt[lamda] Sqrt[x] BesselK[1, Sqrt[2] Sqrt[...


7

Aha, finally here comes a chance to share my implmentation of 2nd order absorbing boundary condition (ABC). To understand why the b.c. is implemented in this way, one may refer to Chapter 6 of Lloyd N. Trefethen, Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations but to be honest I myself don't understand the underlying ...


7

@MichaelE2 gave the idea to use a shooting method, because NDSolve is only able to handle initial value problems involving DiracDelta First solve the problem with a parametric slope f'[0]==fs0 F = ParametricNDSolveValue[{f''[x] + x f'[x] == x^2*DiracDelta[x - 1], f[0] == 0, f'[0] == fs0}, f , {x, 0, 10}, fs0] Now choose fs0 to fullfill the second boundary ...


7

ClearAll[x, y, t, f]; f = 10*Sin[Pi*x]^10*Sin[Pi*y]^10*Sin[Pi*t]^10; pde = D[u[x, y, t], t] == D[u[x, y, t], {x, 2}] + D[u[x, y, t], {y, 2}] + f; ic = u[x, y, 0] == 0; bc = {u[0, y, t] == 0, u[1, y, t] == 0, u[x, 0, t] == 0, u[x, 1, t] == 0}; sol = NDSolve[{pde, ic, bc}, u, {x, 0, 1}, {y, 0, 1}, {t, 0, 4}]; Manipulate[ Plot3D[Evaluate[u[x, y, t0] /. sol],...


7

As mentioned user21 for setting traction boundary conditions we need to know normal and tangent to the surface which is subjected to loading. On the inner surface ($r=a$) for the unit tangent vector we have $\vec{\tau}=(-y/a,x/a)$ that corresponds to anticlockwise direction. On the inner surface the normal component of traction vector is zero, whereas shear ...


7

Based on the suggestion by Daniel Lichtblau, the functional differential equation in the question can be converted to a delay differential equation as follows. With rult = t -> 2^-s; 2 t becomes 2^-(s-1), as desired. To proceed, ruls = Simplify[Solve[Equal @@ rult, s] /. C[1] -> 0, t > 0][[1, 1]] (* s -> -Log[t]/Log[2] *) and the range of ...


7

This is a well-known problem of a rigid domain wall called after Zhirnov, who solved it first. The main point here is that the domain wall solution (that you are, in fact, looking for) corresponds to a separatrix. That is, it is a special trajectory. To select it, one needs to take the corresponding boundary conditions. It brings up the desired simple ...


Only top voted, non community-wiki answers of a minimum length are eligible