MATLAB: Eval returns NaN when the number “r” gets too large

eval() returns nan when my matrix gets really bigthe evil eval

Good morning, i am doing a numerical integration and I need to evaluate the equation below for r = 0.538469810105683 which is a large number but necessary for accurate solution. When I verify the result turns out to be very large and it just appears NaN. What should I do, please any help will be very supporting. Thanks in advance

g(r)=(12*(((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 – (12*r + 12)^3/144115188075855872 – 1688849860263933/281474976710656)^2 + 1)^(1/2)*((15554*r^6)/625 – (19998*r^7)/625 + (140007*r^5)/5000 – (20001*r^4)/1000 – 5*r^3 + 3*r^2 + (417*r)/2500 – 139/2500) – (12*((27*r)/1125899906842624 – (9*(12*r + 12)^2)/36028797018963968 + 39/1125899906842624)^2*((1111*r^9)/2500 – (1111*r^8)/2500 – (6667*r^7)/10000 + (6667*r^6)/10000 + r^5/4 – r^4/4 – (139*r^3)/5000 + (139*r^2)/5000))/(((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 – (12*r + 12)^3/144115188075855872 – 1688849860263933/281474976710656)^2 + 1)^(1/2) – (24*((27*r)/1125899906842624 – (9*(12*r + 12)^2)/36028797018963968 + 39/1125899906842624)*((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 – (12*r + 12)^3/144115188075855872 – 1688849860263933/281474976710656)*((9999*r^8)/2500 – (2222*r^7)/625 – (46669*r^6)/10000 + (20001*r^5)/5000 + (5*r^4)/4 – r^3 – (417*r^2)/5000 + (139*r)/2500))/(((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 – (12*r + 12)^3/144115188075855872 – 1688849860263933/281474976710656)^2 + 1)^(1/2) + (12*((27*r)/1125899906842624 – (9*(12*r + 12)^2)/36028797018963968 + 39/1125899906842624)^2*((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 – (12*r + 12)^3/144115188075855872 – 1688849860263933/281474976710656)^2*((1111*r^9)/2500 – (1111*r^8)/2500 – (6667*r^7)/10000 + (6667*r^6)/10000 + r^5/4 – r^4/4 – (139*r^3)/5000 + (139*r^2)/5000))/(((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 – (12*r + 12)^3/144115188075855872 – 1688849860263933/281474976710656)^2 + 1)^(3/2) + (12*((81*r)/1125899906842624 + 27/562949953421312)*((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 – (12*r + 12)^3/144115188075855872 – 1688849860263933/281474976710656)*((1111*r^9)/2500 – (1111*r^8)/2500 – (6667*r^7)/10000 + (6667*r^6)/10000 + r^5/4 – r^4/4 – (139*r^3)/5000 + (139*r^2)/5000))/(((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 – (12*r + 12)^3/144115188075855872 – 1688849860263933/281474976710656)^2 + 1)^(1/2))^2/(1728*(((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 – (12*r + 12)^3/144115188075855872 – 1688849860263933/281474976710656)^2 + 1)^(3/2))

Best Answer

I am not certain what the problem is.

To evaluate the ‘g’ function at the chosen value of ‘r’, try this:

syms r 
g(r)=(12*(((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 - (12*r + 12)^3/144115188075855872 - 1688849860263933/281474976710656)^2 + 1)^(1/2)*((15554*r^6)/625 - (19998*r^7)/625 + (140007*r^5)/5000 - (20001*r^4)/1000 - 5*r^3 + 3*r^2 + (417*r)/2500 - 139/2500) - (12*((27*r)/1125899906842624 - (9*(12*r + 12)^2)/36028797018963968 + 39/1125899906842624)^2*((1111*r^9)/2500 - (1111*r^8)/2500 - (6667*r^7)/10000 + (6667*r^6)/10000 + r^5/4 - r^4/4 - (139*r^3)/5000 + (139*r^2)/5000))/(((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 - (12*r + 12)^3/144115188075855872 - 1688849860263933/281474976710656)^2 + 1)^(1/2) - (24*((27*r)/1125899906842624 - (9*(12*r + 12)^2)/36028797018963968 + 39/1125899906842624)*((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 - (12*r + 12)^3/144115188075855872 - 1688849860263933/281474976710656)*((9999*r^8)/2500 - (2222*r^7)/625 - (46669*r^6)/10000 + (20001*r^5)/5000 + (5*r^4)/4 - r^3 - (417*r^2)/5000 + (139*r)/2500))/(((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 - (12*r + 12)^3/144115188075855872 - 1688849860263933/281474976710656)^2 + 1)^(1/2) + (12*((27*r)/1125899906842624 - (9*(12*r + 12)^2)/36028797018963968 + 39/1125899906842624)^2*((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 - (12*r + 12)^3/144115188075855872 - 1688849860263933/281474976710656)^2*((1111*r^9)/2500 - (1111*r^8)/2500 - (6667*r^7)/10000 + (6667*r^6)/10000 + r^5/4 - r^4/4 - (139*r^3)/5000 + (139*r^2)/5000))/(((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 - (12*r + 12)^3/144115188075855872 - 1688849860263933/281474976710656)^2 + 1)^(3/2) + (12*((81*r)/1125899906842624 + 27/562949953421312)*((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 - (12*r + 12)^3/144115188075855872 - 1688849860263933/281474976710656)*((1111*r^9)/2500 - (1111*r^8)/2500 - (6667*r^7)/10000 + (6667*r^6)/10000 + r^5/4 - r^4/4 - (139*r^3)/5000 + (139*r^2)/5000))/(((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 - (12*r + 12)^3/144115188075855872 - 1688849860263933/281474976710656)^2 + 1)^(1/2))^2/(1728*(((3*r)/281474976710656 + (3*(12*r + 12)^2)/36028797018963968 - (12*r + 12)^3/144115188075855872 - 1688849860263933/281474976710656)^2 + 1)^(3/2));
r = 0.538469810105683;
Result = vpa(g(r))

producing:

Result =
0.00014765654531859441862062868960765

If you want to use that number outside of the Symbolic Math Toolbox, use the double function:

Result = double(g(r))

producing:

Result =
     1.476565453185944e-04

If you are doing something else, please be specific, and preferably post the code that is causing problems.

Related Solutions

MATLAB: How to applay y”’+y’=0 , y(0)=y0=0, y'(0)=y01=1, y”(0)=y011=-2 in the following system where fn=y”’ i got this error Error using mupadengine/feval (line 157) MuPAD error: Error: Cannot differentiate the equation. [numeric::fsolve] Error in sym

Change

syms y

syms y(t)

MATLAB: When I run the program inaccurate results appear What is the problem

The best approach to this is likely not the vpasolve function, since it takes forever.

I have no idea of these are the ‘accurate’ results you want (I have no idea what that means), however you will get them much more quickly:

syms t y1 y2 y3 y4  yn11 yn1 yr yr11 yr111   yr v y(v) x
format longe
h=0.1;
r=1.7686999343761197980595002942562;   %1694/911=1.859495060373216e+00
x=0;
y0=0;
y01=0.25;
y011=0;
fn=-y01+y0.*y01.*y011;
fnr=-yr11+yr.*yr11.*yr111;
fn1=-yn1+y1.*yn1.*yn11;
fn2=-y3+y2.*y3.*y4;
eq1 =-y2-(r - 2)/r*y0+ 2/(r*(r - 1))*yr+ (2*r - 4)/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 - 23*r + 6))/(120*r)*fn+ (h.^3.*(- r^2 + 2*r + 3))/(60*r)*fnr -(h.^3.*(- r^3 + 4*r^2 + 9*r - 26))/60*fn1+ (h.^3.*(- r^3 + 2*r + 1))/120*fn2;
eq1 = simplifyFraction(eq1);
eq1d = vpa(eq1, 5)
eq2= -h.*y3 -(r - 3)/r*y0+ 3/(r*(r - 1))*yr+ (r - 4)/(r - 1)*y1 -(h.^3.*(3*r^4 - 24*r^3 + 66*r^2 - 67*r + 12))/(240*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r - 4))/(40*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- 3*r^4 + 15*r^3 + 15*r^2 - 137*r + 116))/(120*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 + 3*r - 12))/(80*(r - 2))*fn2;
eq2 = simplifyFraction(eq2);
eq2d = vpa(eq2, 5)
eq3=-h.^2.*y4+2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1+ (h.^3.*(- r^4 + 8*r^3 - 22*r^2 + 18*r + 5))/(120*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r + 5))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- r^4 + 5*r^3 + 5*r^2 - 75*r + 77))/(60*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 + 42*r - 81))/(120*(r - 2))*fn2;
eq3 = simplifyFraction(eq3);
eq3d = vpa(eq3, 5)
eq4=-h.*yn1 -(r - 1)/r*y0+ 1/(r*(r - 1))*yr+ (r - 2)/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 - 21*r + 6))/(240*r)*fn+ (h.^3.*(- r^2 + 2*r + 3))/(120*r)*fnr -(h.^3.*(- r^3 + 4*r^2 + 9*r - 8))/120*fn1 -(h.^3.*(r^3 - 2*r + 1))/240*fn2;
eq4 = simplifyFraction(eq4);
eq4d = vpa(eq4, 5)
eq5=-h.^2.*yn11+ 2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 - 28*r + 10))/(120*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r - 10))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- r^4 + 5*r^3 + 5*r^2 - 35*r + 22))/(60*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 - 8*r + 4))/(120*(r - 2))*fn2;
eq5 = simplifyFraction(eq5);
eq5d = vpa(eq5, 5)
eq6=-h.*yr11+ (r - 1)/r*y0+ (2*r - 1)/(r*(r - 1))*yr -r/(r - 1)*y1+ (h.^3.*(2*r^4 - 13*r^3 + 28*r^2 - 22*r + 5))/240*fn+ (h.^3.*(4*r^3 - 15*r^2 + 10*r + 5))/(120*(r - 2))*fnr+ (h.^3.*r*(- 2*r^3 + 7*r^2 + 2*r - 3))/120*fn1+ (h.^3.*r*(2*r^4 - 5*r^3 + 2*r^2 + 2*r - 1))/(240*(r - 2))*fn2;
eq6 = simplifyFraction(eq6);
eq6d = vpa(eq6, 5)
eq7=-h.^2.*yr111+ 2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1+ (h.^3.*(4*r^4 - 22*r^3 + 38*r^2 - 22*r + 5))/(120*r)*fn -(h.^3.*(- 14*r^4 + 55*r^3 - 55*r^2 + 5*r + 5))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(4*r^4 - 15*r^3 + 5*r^2 + 5*r - 3))/(60*(r - 1))*fn1+ (h.^3.*(4*r^4 - 8*r^3 + 2*r^2 + 2*r - 1))/(120*(r - 2))*fn2;
eq7 = simplifyFraction(eq7);
eq7d = vpa(eq7, 5)
eq8=-h.*y01 -(r + 1)/r*y0 -1/(r*(r - 1))*yr+ r/(r - 1)*y1+ (h.^3.*(r^3 - 8*r^2 + 22*r - 5))/240*fn+ (h.^3.*(r^3 - 5*r^2 + 5*r + 5))/(120*(r^2 - 3*r + 2))*fnr+ (h.^3.*r*(- r^3 + 5*r^2 + 5*r - 3))/(120*(r - 1))*fn1 -(h.^3.*r*(- r^3 + 2*r^2 + 2*r - 1))/(240*(r - 2))*fn2;
eq8 = simplifyFraction(eq8);
eq8d = vpa(eq8, 5)
eq9 =-h.^2.*y011+ 2/r*y0+ 2/(r*(r - 1))*yr -2/(r - 1)*y1 -(h.^3.*(r^4 - 8*r^3 + 22*r^2 + 22*r - 5))/(120*r)*fn -(h.^3.*(r^4 - 5*r^3 + 5*r^2 + 5*r + 5))/(60*r*(r^2 - 3*r + 2))*fnr -(h.^3.*(- r^4 + 5*r^3 + 5*r^2 + 5*r - 3))/(60*(r - 1))*fn1+ (h.^3.*(- r^4 + 2*r^3 + 2*r^2 + 2*r - 1))/(120*(r - 2))*fn2;
eq9 = simplifyFraction(eq9);
eq9d = vpa(eq9, 5)
eqsfcn = matlabFunction([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9], 'Vars',{[y1,y2,y3,y4,yn1,yn11,yr,yr11,yr111]})
eqsol = fsolve(eqsfcn, rand(1,9))
eqsolc = num2cell(eqsol);
[y1,y2,y3,y4,yn1,yn11,yr,yr11,yr111] = eqsolc{:}

I used the simplifyFraction function to do just that with your equations, then the fsolve function to do the actual calculations. These changes make it significantly more efficient.

The results I got are:

y1 =
     2.495835158216658e-02
y2 =
     4.966725033036157e-02
y3 =
     2.450145921764798e-01
y4 =
    -4.970810349993776e-02
yn1 =
     2.487509121241510e-01
yn11 =
    -2.496352514349295e-02
yr =
     4.398727143758653e-02
yr11 =
     2.460985502335340e-01
yr111 =
    -4.401564490668385e-02

You must be certain your equations are written as you want them to be in order to get ‘accurate’ values. We have no control over that.

Experiment to get the results you want.

Best Answer

Related Solutions

MATLAB: How to applay y”’+y’=0 , y(0)=y0=0, y'(0)=y01=1, y”(0)=y011=-2 in the following system where fn=y”’ i got this error Error using mupadengine/feval (line 157) MuPAD error: Error: Cannot differentiate the equation. [numeric::fsolve] Error in sym

MATLAB: When I run the program inaccurate results appear What is the problem

Related Question