0
|
4382 Views
|
4 Replies
|
3 Total Likes
View groups...
Share
GROUPS:

# Issue of DSolve function and 0 ^ 0

Posted 10 years ago
 I was surprised, having faced the following issue in a simple use of DSolve function.I try to find solution of a simple boundary problem for diff. equation of 2nd order:a v''(t) + b v(t) = g(t),v(0)=0, v(1)=1.Coefficients a, b are constants and the left part g(t) is a sum of exponents and linear functions with respect to t. a <0, b >0, so the problem has a unique solution and can be easily calculated by hand.In my case a, b, g(t) calculated on a previous stage.Sapshot of the code and error messages: The code:DSolve[{-0.053399786400854396582413499466002135991456034176000000000000921343149.52287874528033 v1''[t]+0.133390997686009255962976031149875400498398006408000000000002368835549.52287874528033 v1[t]==-0.283018867924528301886792+0.005014252567755304845165779349483873526333219876 E^(-2.81465156744174009869712731637825723649652499385228694169 t)+0.008586511345122480854910534410479478741705706287 E^(-1.732050807568877293527446341505872366942805253810158290235 t)+0.25001517408748706819836422234433174826 E^(-0.80570784117217514973733627569314947705823994352092046887 t)+0.09989713599393614000262596240276808141 E^(0.80570784117217514973733627569314947705823994352092046887 t)-0.028611431245442879573315596710230279738501719103 E^(1.732050807568877293527446341505872366942805253810158290235 t)-0.00248797240353949560222640670452755712872571226127.907709483032015 E^(2.8146515674417400986971273163782572364965249938522869416940.50024606432001 t)+0.2402990388038447846208607475970096119627.907212282747974 t, v1 == 0, v1 == 1}, {v1[t]}, t]But Mathematica fails to solve the problem, displaying a number of errors: Det::mindet: Input matrix contains an indeterminate entry. >>FullSimplify::infd: Expression 0.^(-0. t) simplified to Indeterminate. >>FullSimplify::infd: Expression 0.^(0. t) simplified to Indeterminate. >>Det::mindet: Input matrix contains an indeterminate entry. >>What does it mean?
4 Replies
Sort By:
Posted 10 years ago
 First a suggestion: you might use NDSolve instead of DSolve.To follow up on Sean Clarke's comment, here is a way to replace Real numbers by symbolic coefficients.  One issue is that two of the numbers are very close.  They agree up to 40 digits, but I left them as different coefficients, even though it might slow down DSolve a little. (* store IVP in a variable *) ivp = {-0.053399786400854396582413499466002135991456034176000000000000921343149.52287874528033 v1''[t] +       0.133390997686009255962976031149875400498398006408000000000002368835549.52287874528033 v1[t] == -0.283018867924528301886792 +       0.005014252567755304845165779349483873526333219876 E^(-2.81465156744174009869712731637825723649652499385228694169 t) +       0.008586511345122480854910534410479478741705706287 E^(-1.732050807568877293527446341505872366942805253810158290235 t) +       0.25001517408748706819836422234433174826 E^(-0.80570784117217514973733627569314947705823994352092046887 t) +       0.09989713599393614000262596240276808141 E^(0.80570784117217514973733627569314947705823994352092046887 t) -       0.028611431245442879573315596710230279738501719103 E^(1.732050807568877293527446341505872366942805253810158290235 t) -       0.00248797240353949560222640670452755712872571226127.907709483032015 *        E^(2.8146515674417400986971273163782572364965249938522869416940.50024606432001 t) +      0.2402990388038447846208607475970096119627.907212282747974 t,    v1 == 0, v1 == 1};numbers = Union[Cases[ivp, x_Real :> Abs[x], Infinity]]  (* get the numerical coefficients *)(*Out= {0.002487972403539495602226406705,   0.005014252567755304845165779349483873526333219876,   0.008586511345122480854910534410479478741705706287,   0.02861143124544287957331559671023027973850171910,   0.05339978640085439658241349946600213599145603417600,   0.09989713599393614000262596240276808141,   0.13339099768600925596297603114987540049839800640800,   0.2402990388038447846208607476,   0.25001517408748706819836422234433174826,  0.283018867924528301886792,   0.80570784117217514973733627569314947705823994352092046887,   1.732050807568877293527446341505872366942805253810158290235,   2.8146515674417400986971273163782572364965,   2.8146515674417400986971273163782572364965249938522869417}*)coeffRules =       (* create replacement rules to replace number by coeff[n] *)  Thread /@ {#, Distribute[-#, Rule, Times]} &[numbers -> Array[coeff, Length@numbers]] // Flatten;ivp /. coeffRules  (* see if the rules work *)(*Out= {coeff v1[t] -    coeff (v1^\[Prime]\[Prime])[t] == -E^(t coeff) coeff +    E^(-t coeff) coeff + E^(-t coeff) coeff -    E^(t coeff) coeff + E^(t coeff) coeff + t coeff +    E^(-t coeff) coeff - coeff, v1 == 0, v1 == 1}*)solCoeff = DSolve[ivp /. coeffRules, {v1[t]}, t];  (* solve the symbolic IVP *)Reverse /@ coeffRules[[;; Length@numbers]]         (* show that reversing the rules creates the inverse substitution *)(*Out= {coeff -> 0.002487972403539495602226406705,   coeff -> 0.005014252567755304845165779349483873526333219876,   coeff -> 0.008586511345122480854910534410479478741705706287,   coeff -> 0.02861143124544287957331559671023027973850171910,   coeff -> 0.05339978640085439658241349946600213599145603417600,   coeff -> 0.09989713599393614000262596240276808141,   coeff -> 0.13339099768600925596297603114987540049839800640800,   coeff -> 0.2402990388038447846208607476,   coeff -> 0.25001517408748706819836422234433174826,   coeff -> 0.283018867924528301886792,   coeff -> 0.80570784117217514973733627569314947705823994352092046887,   coeff -> 1.732050807568877293527446341505872366942805253810158290235,   coeff -> 2.8146515674417400986971273163782572364965,  coeff -> 2.8146515674417400986971273163782572364965249938522869417}*)sol = solCoeff /. (Reverse /@ coeffRules[[;; Length@numbers]]);  (* substitute the numbers for the coeff[n] *)Plot[v1[t] /. First@sol, {t, 0, 1}] Comparison with Ilian Gachevski's solutionThere's no significant differencesolFSimp = DSolve[FullSimplify@ivp, {v1[t]}, t];diffFSimp = Table[(v1[t] /. First@sol) - (v1[t] /. First@solFSimp), {t, 0, 1, 1/50000}];{Min[diffFSimp], Max[diffFSimp]}(*Out= {0.*10^-13, 0.*10^-12}*)Comparison with NDSolveWe'll construct two solution, a machine-precision one nsol and a one that takes advantage of the high-precision coefficients.  The machine-precision solution is calculated very fast, and the higher precision solution is slower but takes only a fraction of a second. nsol = NDSolve[ivp, {v1[t]}, {t, 0, 1}];  Min[Precision /@ numbers]  (* find the min. precision of the coefficients *) (* Out=   23.4518 *)  nsolprec = NDSolve[ivp, {v1[t]}, {t, 0, 1},   PrecisionGoal -> Min[Precision /@ numbers] - 3,  (* need some slop in the goals *)   AccuracyGoal -> Min[Precision /@ numbers] - 1,    WorkingPrecision -> Min[Precision /@ numbers]];The default AccuracyGoal and PrecisionGoal are one-half of MachinePrecision, or about 10^-8:Plot[Evaluate[(v1[t] /. First@sol) - (v1[t] /. First@nsol)], {t, 0, 1}, PlotRange -> All] AccuracyGoal and PrecisionGoal are not always achievable, but they make the solution agree much better with the symbolic solution:diffPrec =  Table[(v1[t] /. First@sol) - (v1[t] /. First@nsolprec), {t, 0, 1, 1/50000}];{Min[diffPrec], Max[diffPrec]}(*Out= {-5.3365*10^-16, 2.117*10^-17}*)ListPlot[diffPrec, DataRange -> {0, 1}, PlotRange -> All] Posted 10 years ago
 Thank you, friends. Your tips are very valuable. FullSimplify of the right part of equation solved the problem, at least, for the current moment. I have not yet  tested Rationalize, but will keep this workaround in mind. Regards.
Posted 10 years ago
 Using floating point numbers in DSolve is a bit of a risk.  Symbolic functions like DSolve are written for symbolic values, not floating point numbers.  Although these symbolic functions often do just fine with floating point numbers,  they do transformations that are symbolically correct but sometimes numerically unstable. It's best to be careful when mixing symbolic algorithms with finite precision things. You can always use Rationalize to convert floating point numbers into rational numbers which are symbolic to avoid the problem. When solving a differential equation, I usually solve the equation with undefined parameters and then later subsituted the values of the parameters into the solution instead.
Posted 10 years ago
 What happens is that an intermediate computation yields a close to zero number with very low precision, and it turns into Indeterminate, not unlike 0^0. For a workaround, try DSolve[{FullSimplify[eqn], v1 == 0, v1 == 1}, v1[t], t]