22
Phew, this became more complex than I thought because the objective is not quadratic (as I had first expected). Here some preparations for an efficient evaluation of the energy and its first two derivatives:
Module[{d, PP, P, e, energy, L},
d = 3;
PP = Table[Compile`GetElement[P, Compile`GetElement[e, i], j], {i, 1, 2}, {j, 1, d}];
energy = (Sqrt[Dot[...
6
Here's a way based on this answer to monitor the time integration and compare a step with the previous step. Integration stops when two steps do not change much compare to the change in time. I used the error norm used by NDSolve (see this and this tutorial), since norm < 1 is the bound on the local error of an integration step. One can substitute one'...
5
Similar to @Nasser's answer one could define an error function which compares the solutions u[t,x] for different times t:
Du = 1;
alpha = 4;
T = 3;
pde = D[uind[t, x], t] ==Du*D[uind[t, x], x, x] - D[uind[t, x], x] - alpha;
bc = {uind[t, 0] == 3, (D[uind[t, x], x] /. x -> 1) == 0};
ic = uind[0, x] == 3;
U = NDSolveValue[{pde, ic, bc}, uind, {t, 0, T}, {x, ...
5
btw, you showed no ic, so I made one up.
One way to show when equilibrium is reached, since the right side is free (from BC), is to to see when the solution at the right side stops changing. Like this
ClearAll[uind, x, t];
Du = 1;
alpha = 4;
T = 3;
pde = D[uind[t, x], t] == Du*D[uind[t, x], x, x] - D[uind[t, x], x] - alpha;
bc = {uind[t, 0] == 3, (D[uind[t, ...
5
Not foolproof, but RootApproximant works for this case:
sol = Solve[Log[2*Sqrt[2+Sqrt[3]],x^2-2 x-2]==Log[2+Sqrt[3],x^2-2 x-3],x,Reals] //Values;
ToRadicals@RootApproximant[sol]
{{1 - Sqrt[11 + 4 Sqrt[3]]}, {1 + Sqrt[11 + 4 Sqrt[3]]}}
3
It looks like a bug. Fit tries to create a numerical function with
Experimental`CreateNumericalFunction[{{x, y}},
{cheb[0, x] cheb[0, y], cheb[0, x] cheb[1, y], cheb[0, y] cheb[1, x],
cheb[1, x] cheb[1, y]}, {4}]
(* Experimental`NumericalFunction::dimsl error *)
But it should call it without the extra {} in the first argument:
Experimental`...
3
Just add WorkingPrecision->5 in your Plotcommand
Plot[{(E^(2 h S x + S x^2 -
2 h S x^2) theta (Erfi[((-1 + h) Sqrt[S])/Sqrt[-1 + 2 h]] -
Erfi[(Sqrt[S] (-h - x + 2 h x))/
Sqrt[-1 + 2 h]]))/((Erfi[((-1 + h) Sqrt[S])/
Sqrt[-1 + 2 h]] + Erfi[(h Sqrt[S])/Sqrt[-1 + 2 h]]) x (1 -
x))}, {t, 100, 300}, PlotRange -&...
2
Clear["Global`*"]
$Version
(* "12.3.1 for Mac OS X x86 (64-bit) (June 19, 2021)" *)
Use exact values for known constants and add the option WorkingPrecision to the plot.
theta = 1/10;
ep = 150;
w2 = 1/10;
c = 1/2;
rho[t_] := 1 + c Sin[w2 t];
h = 1/4;
S = ep (1 + c Sin[w2 t]);
x = 1/5;
Plot[{(E^(2 h S x + S x^2 -
2 h S x^2) ...
2
You might consider a single function for calculating the entropy values:
e[y_?NumericQ, d_] := Module[{skd, dmin, dmax, g},
skd = SmoothKernelDistribution[d, y];
{dmin, dmax} = MinMax[d] + 3.9 y {-1, 1};
g = PDF[skd, a][[1, 1, 1]];
NIntegrate[-g Log[g], {a, dmin, dmax}]]
Rather than integrating from $-\infty$ to $\infty$, I've used 3.9 times the ...
2
Since you are using NDSolve to find an approximate solution to your pde, maybe an approximate Dirac Delta will work for you. Since a true Dirac Delta is an infinitely high, zero width spike with an area of one, I use:
dd[x_] = aa E^(-100000 x^2)
Solve for aa to make the area 1.
Integrate[dd[x], {x, -Infinity, Infinity}] == 1
aa = aa /. Solve[%, aa][[1]]
(*...
2
I wonder if the issue is just that of lack of numerical precision and/or there are multiple solutions that give the same minimum. Consider the pieces of your code as error (the function being minimized), constraints, and vars.
Using NMinimize:
results = Table[10^i NMinimize[{error/10^i, constraints}, vars][[1]], {i, 6, 10}]
[![Unbounded error message][1]][1]...
2
Using
QuadraticOptimization[error, constraints, vars, Method -> "COIN",
Tolerance -> 0.000001]
(possibly tweaking the Tolerance) instead gives good answers.
1
Here is another method that is more general in one way and less in another. It uses the values of the time derivative as the measure of when equilibrium has been reached. I approximated the derivative in my other answer, but since the derivative is computed at each step by NDSolve anyway, it makes sense to try to use it. OTOH, it uses the MonitorMethod (...
1
Here is a different approach that may raise some issues with the Question itself. As I noted in comments above, the functions listed in the question are quite noisy, and FindMinimum does not converge in any reasonable amount of time. Consequently, as suggested by the OP in a comment above, JimB in his answer used an alternative function,
e[y_?NumericQ, d_] ...
1
Plotting f suggests that it has very many, if not an infinite number of zeros.
Plot[f, {x, 0, Pi}, PlotPoints -> 10000, ImageSize -> Large,
AxesLabel -> {x, "f"}, LabelStyle -> {15, Bold, Black}]
Plot[f, {x, 3, Pi}, PlotPoints -> 100000, ImageSize -> Large,
AxesLabel -> {x, "f"}, LabelStyle -> {15, Bold, ...
1
Combining the FindRoot usage suggested by @DerekH with manipulate:
Manipulate[
sol = FindRoot[{eq1, eq2}, {{x, xcoord}, {y, ycoord}}];
Show[{ContourPlot[Evaluate@{eq1, eq2}, {x, -10, 10}, {y, -10, 10},
PlotLabel -> Style[sol, Red, Bold], GridLines -> Automatic],
Graphics[{Red, PointSize[Medium], Point[{x, y} /. sol]}]
}],
{{xcoord, 0}, 0, ...
1
Restricting y too, we succeed:
NSolve[{eq1, eq2, x >= 0.01, x <= 0.5, y >= -5, y <= 5}, {x, y}, Reals]
or
NSolve[{eq1, eq2, x \[Element] Interval[{0.01, 0.5}], y >= -5, y <= 5}, {x, y}, Reals]
{{x -> 0.0117594, y -> 0.418582}}
1
Since the real part of the integrand is symmetric about t1==t2, simply double the integral and integrate t2 from zero to t1 with no problems.The imaginary part is antisymmetric and disappears. Therefore work only with Cos
By the way, this is no improper integral, as limit t2->t1 is a finite number. Limit[\[Delta]Et\[Delta]E0[t1, t2, 1, 1]*E^(I*1*(t1 - t2))...
Only top voted, non community-wiki answers of a minimum length are eligible