## NDSolve::ndsz: step size is effectively zero; singularity or stiff system suspected + other warnings for system of differential equations

I’m trying to solve a set of differential equations numerically to get a 3D plot, but I am getting multiple different warnings and errors. First of all, here is the code:

ClearAll["Global*"] tdot[\[Lambda]_, q_] := Simplify[1/(1 - 1/r)] /. r -> r[\[Lambda]] rdot[\[Lambda]_, q_] := Simplify[(1/r)*Sqrt[r^2 - q^2*(1 - 1/r)]] /. r -> r[\[Lambda]] \[Phi]dot[\[Lambda]_, q_] := Simplify[-(q/r^2)] /. r -> r[\[Lambda]] tasym[\[Lambda]_, q_] := \[Lambda] + Log[\[Lambda]] - 1/\[Lambda] + (q^2 - 2)/(4*\[Lambda]^2) + (3*q^2 - 4)/(12*\[Lambda]^3) + (-3*q^4 + 8*q^2 - 8)/(32*\[Lambda]^4) +     (-9*q^4 + 20*q^2 - 16)/(80*\[Lambda]^5) rasym[\[Lambda]_, q_] := \[Lambda] + q^2/(2*\[Lambda]) - q^2/(4*\[Lambda]^2) - q^4/(8*\[Lambda]^3) + (3*q^4)/(16*\[Lambda]^4) + (q^6/16 - q^4/20)/\[Lambda]^5 \[Phi]asym[\[Lambda]_, q_] := q/\[Lambda] - q^3/(3*\[Lambda]^3) + q^3/(8*\[Lambda]^4) + q^5/(5*\[Lambda]^5) \[Lambda]min = 0;  \[Lambda]max = 20;  \[Lambda]inf = 999;  qlist = Array[N[#1] & , 100, {-20, 20}];  maxder = 999999;  eq[q_] := {t[\[Lambda]], r[\[Lambda]], \[Phi][\[Lambda]]} /. NDSolve[{Derivative[1][t][\[Lambda]] == tdot[\[Lambda], q], Derivative[1][r][\[Lambda]] == rdot[\[Lambda], q],        Derivative[1][\[Phi]][\[Lambda]] == \[Phi]dot[\[Lambda], q], t[\[Lambda]inf] == tasym[\[Lambda]inf, q], r[\[Lambda]inf] == rasym[\[Lambda]inf, q],        \[Phi][\[Lambda]inf] == \[Phi]asym[\[Lambda]inf, q], WhenEvent[{Abs[Derivative[1][t][\[Lambda]]] > maxder ||           Abs[Derivative[1][r][\[Lambda]]] > maxder || Abs[Derivative[1][\[Phi]][\[Lambda]]] > maxder}, "StopIntegration"]},       {t[\[Lambda]], r[\[Lambda]], \[Phi][\[Lambda]]}, {\[Lambda], \[Lambda]min, \[Lambda]inf},       {"ExtrapolationHandler" -> {Indeterminate & , "WarningMessage" -> False}}][[1]] eqlist = (eq[#1] & ) /@ qlist;  tlist = eqlist[[All,1]];  rlist = eqlist[[All,2]];  \[Phi]list = eqlist[[All,3]];  surface = MapThread[{#2*Sin[#3], #2*Cos[#3], If[Abs[#3] < Pi, #1, Indeterminate]} & , {tlist, rlist, \[Phi]list}];  Show[ParametricPlot3D[surface, {\[Lambda], \[Lambda]min, \[Lambda]max}, PlotRange -> {{-20, 20}, {-30, 20}, {-30, 20}}], ImageSize -> Large] 

Running this code, with the parameters I chose after playing around with them, gives no warning message, and results in what i’m trying to get, but only half of it. For certain values of the parameters $$\lambda_{inf}$$, $$maxder$$, the warning message in the title of the question appears, however after reading other questions regarding this issue I don’t think this is too much of a problem.

The main problems start when i set $$\lambda_{min}=-30$$ which would give me the other half of the plot. the first warning message that appears is

NDSolve::mxst: Maximum number of 129336 steps reached at the point [Lambda] == -0.53469.

same thing for other values of $$\lambda$$. At first I tried overcoming this by increasing MaxSteps, however this didn’t work for me as Mathematica would just use up all my RAM and my computer would stop working.

To investigate I tried to solve only for the negative values, I set $$\lambda_{min}=-30$$ and $$\lambda_{max}=-5$$, and set NDSolve to solve only between $$\lambda_{min}$$ and $$\lambda_{max}$$. What happens with these settings is the following warning message:

NDSolveProcessSolutions::nodata: No solution data was computed between [Lambda] == -30. and [Lambda] == -5..

Which is weird since I should have a solution, unless I made some dumb mistake. Reading other questions, I saw this could maybe be a case of backslide where since the solution doesn’t adhere to the new standards of NDSolve it doesn’t give any solution. However this is just speculation.

I also tried adding some methods to NDSolve with no results, since I don’t know much about them. What I hope is that by tweaking some NDSolve parameters or using some methods I can manage to get a result, so any suggestion in this direction is very welcome.

## How do I write math equations in my blog?

I have been trying to do this for a long time and none of the methods I have come across seem to work for me. The most common suugestion I have seen is to use MathJax which can be activated by using this script

    <script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script> <script async="" id="MathJax-script" src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script> 

This still doesn’t seem to work for me or maybe I am just doing it the wrong way. I know how to use LATEX but i am not familiar with html.

## solving self consistent equations

I couldn’t attach pic of my self consistent equations in the comment section that I am trying to solve in an earlier message threading

(Elimination of the variables)

I am trying to solve for z for given x = 0.057, y = 0.01 self consistently

y = 1/100;  x = 57/1000;  gE[(n_Integer)?NonNegative] := 1/g[n] - Sqrt[1 + z^2/(Pi*(2*n + 1)*y + (44*x*Pi*g[n])/10)^2];  sce[0] = (2*Pi*1*((Pi*1*(g[0] - 1))/(100*(((Pi*1)/100 + (44*Pi*x)/10)*((Pi*1)/100 + (44*Pi*x*g[0])/10)))))/100 + 6431/10000;  sce[(n_Integer)?Positive] := (2*Pi*1*(((2*n + 1)*Pi*1*(g[n] - 1))/(100*((((2*n + 1)*Pi*1)/100 + (44*Pi*x)/10)*(((2*n + 1)*Pi*1)/100 + (44*Pi*x*g[n])/10)))))/100;  nmax = 3;  nValues = Range[0, nmax];  ee = Flatten[{Total[sce /@ nValues], gE /@ nValues}];  xed0 = GroebnerBasis[ee, {z}, g /@ nValues, MonomialOrder -> EliminationOrder];   NSolve[xed0 == 0, z] 

For nmax = 3 it doesn’t converge and keeps on running and I am aware of @ Daniel Lichtblau comment that here complexity arises as nmax increases.

Question: Is there any other way of solving above self consistent equations?

## System of Differential equations on non uniform grid

I need to solve the following set of ODE’s numerically,

$$\frac{dy}{dx}=Ay+Bz, \ \frac{dz}{dx}=Cy+Dz$$

where the independent variable $$x$$ is an non-uniform spaced array of points and so is A, B, C & D.

What is the best way of achieving this as common NDsolve specifies the independent variable in $$\{x,x_{min},x_{max}\}$$ format.

## Simplifying equations in python using wolframclient

I’ve this mathematical expression and want it to Simplify it in python using wolframclient wrapper.

1/700*((-1316*w + b - 0.0)^2 + (-1116*w + b - 1.0)^2) 

I’ve tried this method:

In[1]: mathematica_session.evaluate(wlexpr('1/700*((-1316*w+b-0.0)^2+(-1116*w+b-1.0)^2)')) Out[1]: Times[Rational[1, 700], Plus[Power[Plus[0.0, Globalb, Times[-1316, Globalw]], 2], Power[Plus[-1.0, Globalb, Times[-1116, Globalw]], 2]]] 

And then this method:

In[2]: mathematica_session.evaluate(wl.Simplify('1/700*((-1316*w+b-0.0)^2+(-1116*w+b-1.0)^2)')) Out[2]: '1/700*((-1316*w+b-0.0)^2+(-1116*w+b-1.0)^2)' 

And finally this is the result when I tried this on directly Mathematica:

In[1]: Simplify[1/700*((-1316*w + b - 0.0)^2 + (-1116*w + b - 1.0)^2)] Out[1]: 0.00285714 (0.5 + b^2 + b (-1. - 2432. w) + 1116. w + 1.48866*10^6 w^2) 

And this is the result I want in python. So what I’m missing or doing wrong in python?

## Getting Mathematica to solve a system of two second order nonlinear ordinary differential equations

I tried solving a system of two second order nonlinear ordinary differential equations using the DSolve command. First, I tried like this:

eqns = {A''[x] == 2/B[x]*A'[x]*B'[x],     B''[x] + 1/B[x]*(A'[x])^2 - 1/B[x]*(B'[x])^2 == 0}; sol = DSolve[eqns, {A, B}, x] 

However, as Mathematica didn’t (couldn’t?) solve this, I transformed it into a system of four first order equations:

eqns = {c'[x] == 2/B[x]*c[x]*d[x],     d'[x] + 1/B[x]*(c[x])^2 - 1/B[x]*(d[x])^2 == 0, c[x] == A'[x],     d[x] == B'[x], c[0] == 1, d[0] == 1, A[0] == 1, B[0] == 1}; sol = DSolve[eqns, {A, B, c, d}, x] 

This still doesn’t work. Weirdly enough, I don’t even get an error message.

I only included the boundary conditions thinking that they may be helpful, but they aren’t part of my original problem.

Your help would be greatly appreciated:)

## Solve Matrix equations with Cross Product: weird system of equations

I would like to find the values {Pfx, Pfy, Pfz} that satisfy the equation A X B = C . The code of everything is at the end, I want to ilustrate with images what i think of first:

A is this:

B is this: {0,0,0}

And C is this: {100,500,200}

The image of the complete code is this one:

My variable here are {Pfx,Pfy,Pfz} and the {i,j,k} are the unit vectors. The way that the matrix is shown in the picture corresponds with a trick used to solve in papper this type of matrix.

The solution would give me the value of Pfx in the "x" coordinate (and this would be expressed by Pfx being multiplied with the i vector). And the same mechanism applies with Pfy related with j Vector ; and Pfz related with k.

The problem here comes with the fact that the I cant find the values {Pfx, Pfy, Pfz} that satisfy the equation A X B = C. I am not sure if the problem lies in the "LinearSolve" comand, in the use of the CrossProduct or in the use of the versor {i, j, k} inside the matrix.

Any kind of help in this regard will be extreamly useful, thanks in advance!!

Code:

i := {1, 0, 0}  j := {0, 1, 0}  k := {0, 0, 1}  LinearSolve[({ {i, j, k},{1, 2, 3},{Pfx, Pfy, Pfz}})\[Cross]({{0},{0},{0}}) == ( {{100},{500},{200}} )]    

## Real solutions of third and fourth degree equations

A few hours ago I "discovered" that if a third or fourth degree equation has distinct real solutions, it’s possible to calculate them avoiding complex numbers.

In particular, we have:

poly = (x - 1)(x - 2)(x - 3); {c, b, a} = CoefficientList[poly, x][[1 ;; 3]]; d = a^2/3 - b; e = 2 a^3/27 - a b/3 + c; f = ArcCos[-3 Sqrt[3] e/(2 d Sqrt[d])]; N[{-a/3 + 2 Sqrt[3 d] Cos[(f - 4 Pi)/3]/3,    -a/3 + 2 Sqrt[3 d] Cos[(f - 2 Pi)/3]/3,    -a/3 + 2 Sqrt[3 d] Cos[(f - 0 Pi)/3]/3}] 

{1., 2., 3.}

poly = (x - 1)(x - 2)(x - 3)(x - 4); {d, c, b, a} = CoefficientList[poly, x][[1 ;; 4]]; e = 3 a^2/4 - 2 b; f = 2 c - a b + a^3/4; g = b^2 + 12 d - 3 a c; h = 27 a^2 d - 9 a b c + 2 b^3 - 72 b d + 27 c^2; i = ArcCos[h/(2 g Sqrt[g])]; j = (e + 2 Sqrt[g] Cos[i/3])/3; N[{-a/4 - (Sqrt[j] + Sqrt[e - j + f/Sqrt[j]])/2,    -a/4 - (Sqrt[j] - Sqrt[e - j + f/Sqrt[j]])/2,    -a/4 + (Sqrt[j] - Sqrt[e - j - f/Sqrt[j]])/2,    -a/4 + (Sqrt[j] + Sqrt[e - j - f/Sqrt[j]])/2}] 

{1., 2., 3., 4.}

After a few moments of joy, however, I noticed that, for example, if I write (x - 10^-15) instead of (x - 1), I get respectively 6.66134*10^-16 and 8.88178*10^-16 instead of 10^-15 as the solution.

I intuitively believe that it’s a numerical problem and that the main cause is ArcCos, but I’m not too sure and also I’m not able to establish if something can be done to solve the issue, or if I have to give up and I must necessarily rely on it to the good Newton method.

Thank you!

## Fit data to a model of two-variable differential equations

I am trying to fit some data to the solution of a differential equation. I am using "ParametricNDSolve" to create a model which should fit the data using "FindFit".

Now, the code I propose is very similar to the examples I have seen to solve similar problems. However, the equation of this case is a little more complicated that the examples of Mathematica manuals (because is a function of two variables, and the initials conditions involves two numeric integrals).

The code keeps running for hours, and it does not appear to solve the problems at all:

data = {{21/100, 0.260882276365091}, {6/25, 0.330910580727009}, {27/ 100, 0.271601829283246}, {3/10, 0.294039749043066}, {33/100,  0.373363616981994}, {9/25, 0.467495450916}, {39/100,  0.50972503848512}, {21/50, 0.639300915026114}, {9/20,  0.679672806174314}, {12/25, 0.693446859556429}, {51/100,  0.697043731283207}, {27/50, 0.736147748563448}, {57/100,  0.706580484286456}, {3/5, 0.719128578284264}, {63/100,  0.762276173639189}, {33/50, 0.788753648395851}, {69/100,  0.802107836002146}, {18/25, 0.803992011056852}, {3/4,  0.796349018322968}, {39/50, 0.795845629380594}, {81/100,  0.808805024532776}, {21/25, 0.812210415525446}, {87/100,  0.811677706096654}, {9/10, 0.796614634731033}, {93/100,  0.802573366818874}, {24/25, 0.801851431134633}, {99/100,  0.799762250401537}, {51/50, 0.803234757269179}, {21/20,  0.812346679044962}, {27/25, 0.804974524275207}, {111/100,  0.813054479870107}, {57/50, 0.802494974488922}, {117/100,  0.811124462792484}, {6/5, 0.810259975308694}, {123/100,  0.810910967136978}, {63/50, 0.815891670404585}, {33/25,  0.803843274767398}, {69/50, 0.814778081357444}, {36/25,  0.809686831953762}, {3/2, 0.829532613157218}, {9/5,  0.815351190536115}, {21/10, 0.828143435644695}, {12/5,  0.791269047273677}, {27/10, 0.821321230878138}, {3,  0.81153212440728}, {33/10, 0.821798500231733}, {18/5,  0.819594302125574}, {39/10, 0.828629593518031}};  model = ParametricNDSolveValue[{ D[v[x, t], t] == d1*D[v[x, t], x, x], v[-5, t] == 0, v[10, t] == flux, v[x, 0] == flux*NIntegrate[Cos[y]^n1, {y, ArcTan[-b1*(x + a1)], Pi/2}]/NIntegrate[Cos[y]^n1, {y, -Pi/2, Pi/2}]}, v, {x, -5, 10}, {t, 0, 2}, {d1, flux, n1, b1, a1}]  fit = FindFit[data,  (c*model[d1, flux, n1, b1, a1][x, 2])/((c - 1)*model[d1, flux, n1, b1, a1][x, 2] + 1), {{d1, 0.001}, {flux, 0.57}, {c, 3.2}, {n1, 2}, {b1,1}, {a1, -0.4}}, x] `

I guess that this code is computationally inefficient, but I need help to understand why.

Thank you all very much

Ps: english is obviously not my first lenguage, so I apologize for the grammar mistakes.

## Pass equations into shaders to define graphics – HLSL or other shaders

Is it possible with HLSL (or other popular shader languages) to pass instead of an image, an equation that would define the pixel color / position output by the shader? This would allow for more dynamic drawing of smooth shapes that could be sampled with infinite detail. They would also by able to change and morph in real time.

Not sure how this is possible without passing a string with the eq to the shader and having it parse it. Seems like that might lose all the performance gained by putting it onto the GPU.