10
Stiffness and singularities
Around 1960, things became completely different and everyone became aware that the world was full of stiff problems. (G. Dahlquist in Aiken 1985)
Shampine's "flame" model $y'=y^2-y^3$, $y(0)=10^{-3}$, $0 \le x \le 2/y(0)$, integrated with a relative precision goal of $10^{-4}$. The red dots show the steps and ...
8
The issue here is that AdjacencyMatrix returns a non-real valued matrix
g = GridGraph[{20, 20}];
amg = AdjacencyMatrix[A1];
meaning the elements are all integers:
In[]:= Map[Head,Normal@amg,{2}]//Flatten//Union
Out[]= {Integer}
and the kernel is trying to find an answer in terms of integers, rational numbers, and roots. Numericizing before computation ...
6
If you just want the largest eigenvalue, the Arnoldi method is much faster than calculating all eigenvalues (and associated eigenvectors) and picking the largest:
Eigensystem[A // N, 1,
Method -> {"Arnoldi", "Criteria" -> "RealPart"}] // AbsoluteTiming
(* {0.007187, {{3.95532}, {{-0.00211558, -0.0041839, ......
5
Changing the numerical integration strategy from the default "GlobalAdaptive" to "LocalAdaptive" gives the expected behavior (and results much closer to the expected ones):
NIntegrate[Exp[-(x - y)^2], {x, -100, 100}, {y, -100, 100}, Method -> "LocalAdaptive"]
(* 352.491 *)
NIntegrate[Exp[-(x - y)^2], {x, -200, 200}, {y, -200,...
5
Try first to simplify, then to substitute:
Simplify[func, {\[Rho] > 0, p > 0, \[Sigma] > 0}] /. {\[Sigma] -> 3./
100, p -> 0.26, \[Rho] -> 0.8}
(* E^(-33.3333 (-0.8 + (x^2 + y^2 +
z^2)^0.26)^2) + E^(-33.3333 (0.8 + (x^2 + y^2 + z^2)^0.26)^2) *)
Have fun!
5
Floats cause that. The Rationalize command helps.
FullSimplify[Rationalize[E^(-33.3333 (0.8 - (x^2 + y^2 + z^2)^0.26)^2) +
E^(-33.3333 (0.8 + (x^2 + y^2 + z^2)^0.26)^2), 10^-40]]
E^(-((333333 (4 - 5 (x^2 + y^2 + z^2)^(13/50))^2)/250000)) + E^(-(( 333333 (4/5 + (x^2 + y^2 + z^2)^(13/50))^2)/10000))
Addition. Maybe, the result of
FullSimplify[E^Expand@...
4
In this case,
?NumericQ condition is not necessary in function definitions, but
?NumericQ is needed for NIntegrate (ref. the documentation for message), so the integral should be split off, here to lambda[z_?NumericQ]:
using := for evaluation to happen only when arguments are passed in.
ClearAll[f, g, h, q, lambda];
f[x_] = Sum[(1 - Abs[k]/6)*E^(I*k*(x - (...
3
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]}
Attempting to solve these ODEs with DSolve is unsuccessful, returning the Solve::ifun error message. That this is the same error obtained by the OP is not surprising, because ...
3
One can do everything analytically:
Simplification of the integral
ir = ComplexExpand[Re[(Cos[x] - h - I*Sin[x])/Abs[Cos[x] - h - I*Sin[x]]]] // Simplify
ii = ComplexExpand[Im[(Cos[x] - h - I*Sin[x])/Abs[Cos[x] - h - I*Sin[x]]]] // Simplify
$$\frac{\cos (x)-h}{\sqrt{h^2-2 h \cos (x)+1}}$$
$$-\frac{\sin (x)}{\sqrt{h^2-2 h \cos (x)+1}}$$
Integration of the ...
3
Note this is not a coupled system of ODE's. The two equations even contradict each other.
You could write f as a function of one variable, then you get the following equations:
f[z_] = -z^2;
eq = {x'[t] + y'[t] == f'[x[t] + y[t]], D[f'[x[t] + y[t]], t] + f'[x[t] + y[t]]}
This is:
The first equations claims: x'[t] + y'[t]== -2(x[t]+y[t]), whereas the ...
3
Two insightful responses were left in comments which I think are worth summarizing here as answers:
As J.M. mentioned, Equal has a tolerance within which the two arguments are considered equal. This is controlled by Internal`$EqualTolerance (see e.g. this question, among others).
Regarding the second question, Michael mentioned that the machine real ...
2
Your function is numerically ill defined around zero. But you may replace the function by a limit. And you need to increase the working precision
f[x0_] := Limit[(1 - Exp[-4*xx^6])/(1 - (1 - xx^6)^(4)) - 1, xx -> x0]
LogLinearPlot[{f[x]}, {x, 10^-6, 1}, PlotRange -> {-.1, .1}, WorkingPrecision -> 20]
2
Is this what you want?
I interpreted: "I am trying to find the maximum value h[ts]
such that my function h is bigger than h[ts] on half of the interval [0,2Ï€] and smaller than h[ts] on the other half, i.e.:"
f[x_] = Sum[(1 - Abs[k]/6)*E^(I*k*(x - (Pi/3))), {k, -5, 5}];
g[x_] = Sum[(1 - Abs[k]/6)*E^(I*k*(x - (5*Pi/3))), {k, -5, 5}];
h[t_] = ...
2
Update
I was complicating matters. No need to split:
ClearAll[f, int, g];
f[x_, h_] := Sign[Cos[x] - h - I Sin[x]]
int[h_] := NIntegrate[f[x, h], {x, 0, 2 Pi}]
g[h_?NumericQ] := Re[2/(1 + 1/(2 Pi) int[h])]
Needs["NumericalCalculus`"]
Plot[ND[g[h], {h, 2}, h0], {h0, 0, 2}]
Taking $69$sec with no message:
Oscillating is perhaps due to numerical ...
2
A method to proceede if there is no analytical solution for the integral and to get very good results.
Generate an interpolating function of h you can easily differentiate where you don't need further numerical integrations when you differentiate.
Therefore differentate the integrand for h, integrate for x, that means i interchanged integration:
(* int[h_] =...
1
Knowing I'm a little bit late, here a straightforward numerical solution.
First define the integral
int = Function[h,NIntegrate[(-h + Cos[x] - I Sin[x])/Sqrt[1 + h^2 - 2 h Cos[x]], {x, 0, 2 Pi}]]
Second define the second derivative
d2 = Function[h,Re@Derivative[2][Function[h, 2/(1 + int[h]/(2 Pi))]][h]]
Plot result
Plot[d2[h], {h, 0, 2}]
1
$Version
(* "12.0.0 for Mac OS X x86 (64-bit) (April 7, 2019)" *)
Clear["Global`*"]
f[x_] = (3 - x^2)^(1/2);
The problem appears to arise because the function is not real throughout the region specified by the constraint 1 < x < 4. Change the constraint to
dom = FunctionDomain[{f[x], 1 < x < 4}, x]
(* 1 < x <= Sqrt[...
1
I assume that x[t] in the definition of L2 is a typo and should read x[r].
With this assumption, you can simply define L2 as:
L2 = r^3 ((V2[x[r], r])*(1 + ((1/2)*(x'[r]^2))) + (1/Pi)*(Pi*(800) -
x[r] + (1/2) Sin[2*x[r]]))
The equations then read:
eq={D[D[L2, x'[r]], r] == D[L2, x[r]], x[0] == 0, x'[0] == 0}
Furthermore, the syntax of NDSolveis ...
Only top voted, non community-wiki answers of a minimum length are eligible
