MATLAB: How to solve a very long non-linear equation with constants

constantssolve equation

I have the following extremely long equation and I need to set it to zero and find the solution in "r"; with "m", "M", and "a" as constants (i.e. the answer would be an equation of r depending on the constants, and not a value):

D2 = 4*m*(2*M – 2*r) – 4*m*r + (2*m^2*(3*r + 2*M*a^2 + a^2*r)*(2*r – 2*M + (M*a)/(2*(M*r)^(1/2)))^2)/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) – (4*m^2*(a^2 + 3)*(a*(M*r)^(1/2) – 2*M*r + r^2)^2)/(r^3*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (6*m^2*(3*r + 2*M*a^2 + a^2*r)*(a*(M*r)^(1/2) – 2*M*r + r^2)^2)/(r^4*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (2*M*m^2*(a^2 – 2*a*(M*r)^(1/2) + r^2)^2)/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (2*m^2*(3*r + 2*M*a^2 + a^2*r)*(2*r – 3*m + (M*a)/(M*r)^(1/2))^2*(a*(M*r)^(1/2) – 2*M*r + r^2)^2)/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^3) – (2*m^2*((M^2*a)/(4*(M*r)^(3/2)) – 2)*(3*r + 2*M*a^2 + a^2*r)*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (4*m^2*(a^2 + 3)*(2*r – 2*M + (M*a)/(2*(M*r)^(1/2)))*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) – (8*m^2*(3*r + 2*M*a^2 + a^2*r)*(2*r – 2*M + (M*a)/(2*(M*r)^(1/2)))*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^3*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) – (4*M*m^2*(2*r – (M*a)/(M*r)^(1/2))*(a^2 – 2*a*(M*r)^(1/2) + r^2))/(r*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (2*M*m^2*(2*r – (M*a)/(M*r)^(1/2))^2*(2*M – r))/(r*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (2*M*m^2*(a^2 – 2*a*(M*r)^(1/2) + r^2)^2*(2*r – 3*m + (M*a)/(M*r)^(1/2)))/(r*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2) + (m^2*((M^2*a)/(2*(M*r)^(3/2)) – 2)*(3*r + 2*M*a^2 + a^2*r)*(a*(M*r)^(1/2) – 2*M*r + r^2)^2)/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2) – (2*m^2*(a^2 + 3)*(2*r – 3*m + (M*a)/(M*r)^(1/2))*(a*(M*r)^(1/2) – 2*M*r + r^2)^2)/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2) + (4*m^2*(3*r + 2*M*a^2 + a^2*r)*(2*r – 3*m + (M*a)/(M*r)^(1/2))*(a*(M*r)^(1/2) – 2*M*r + r^2)^2)/(r^3*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2) + (2*M*m^2*(2*M – r)*(a^2 – 2*a*(M*r)^(1/2) + r^2)^2)/(r^3*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (M*m^2*((M^2*a)/(2*(M*r)^(3/2)) – 2)*(2*M – r)*(a^2 – 2*a*(M*r)^(1/2) + r^2)^2)/(r*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2) – (8*a*m^3*(M*r)^(1/2)*(2*r – (M*a)/(M*r)^(1/2))*(2*r – 2*M + (M*a)/(2*(M*r)^(1/2))))/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) – (4*M*m^2*(2*r – (M*a)/(M*r)^(1/2))*(2*M – r)*(a^2 – 2*a*(M*r)^(1/2) + r^2))/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (16*a*m^3*(M*r)^(1/2)*(2*r – (M*a)/(M*r)^(1/2))*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^3*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (2*M*m^2*(2*M – r)*(a^2 – 2*a*(M*r)^(1/2) + r^2)^2*(2*r – 3*m + (M*a)/(M*r)^(1/2)))/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2) + (4*a*m^3*(M*r)^(1/2)*((M^2*a)/(4*(M*r)^(3/2)) – 2)*(a^2 – 2*a*(M*r)^(1/2) + r^2))/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (16*a*m^3*(M*r)^(1/2)*(a^2 – 2*a*(M*r)^(1/2) + r^2)*(2*r – 2*M + (M*a)/(2*(M*r)^(1/2))))/(r^3*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) – (24*a*m^3*(M*r)^(1/2)*(a^2 – 2*a*(M*r)^(1/2) + r^2)*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^4*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (2*M*m^2*(2*M – r)*(a^2 – 2*a*(M*r)^(1/2) + r^2)^2*(2*r – 3*m + (M*a)/(M*r)^(1/2))^2)/(r*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^3) + (2*M*m^2*((M^2*a)/(2*(M*r)^(3/2)) + 2)*(2*M – r)*(a^2 – 2*a*(M*r)^(1/2) + r^2))/(r*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) – (4*m^2*(3*r + 2*M*a^2 + a^2*r)*(2*r – 2*M + (M*a)/(2*(M*r)^(1/2)))*(2*r – 3*m + (M*a)/(M*r)^(1/2))*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2) – (4*a*m^3*(M*r)^(1/2)*((M^2*a)/(2*(M*r)^(3/2)) + 2)*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (8*a*m^3*(M*r)^(1/2)*(a^2 – 2*a*(M*r)^(1/2) + r^2)*(2*r – 2*M + (M*a)/(2*(M*r)^(1/2)))*(2*r – 3*m + (M*a)/(M*r)^(1/2)))/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2) – (16*a*m^3*(M*r)^(1/2)*(a^2 – 2*a*(M*r)^(1/2) + r^2)*(2*r – 3*m + (M*a)/(M*r)^(1/2))*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^3*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2) + (M^2*a*m^3*(a^2 – 2*a*(M*r)^(1/2) + r^2)*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^2*(M*r)^(3/2)*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) – (4*M*a*m^3*(2*r – (M*a)/(M*r)^(1/2))*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^2*(M*r)^(1/2)*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) – (8*a*m^3*(M*r)^(1/2)*(a^2 – 2*a*(M*r)^(1/2) + r^2)*(2*r – 3*m + (M*a)/(M*r)^(1/2))^2*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^3) – (4*M*m^2*(2*r – (M*a)/(M*r)^(1/2))*(2*M – r)*(a^2 – 2*a*(M*r)^(1/2) + r^2)*(2*r – 3*m + (M*a)/(M*r)^(1/2)))/(r*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2) + (8*a*m^3*(M*r)^(1/2)*(2*r – (M*a)/(M*r)^(1/2))*(2*r – 3*m + (M*a)/(M*r)^(1/2))*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2) – (4*M*a*m^3*(a^2 – 2*a*(M*r)^(1/2) + r^2)*(2*r – 2*M + (M*a)/(2*(M*r)^(1/2))))/(r^2*(M*r)^(1/2)*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) – (4*a*m^3*(M*r)^(1/2)*((M^2*a)/(2*(M*r)^(3/2)) – 2)*(a^2 – 2*a*(M*r)^(1/2) + r^2)*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^2*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2) + (8*M*a*m^3*(a^2 – 2*a*(M*r)^(1/2) + r^2)*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^3*(M*r)^(1/2)*(2*a*(M*r)^(1/2) – 3*m*r + r^2)) + (4*M*a*m^3*(a^2 – 2*a*(M*r)^(1/2) + r^2)*(2*r – 3*m + (M*a)/(M*r)^(1/2))*(a*(M*r)^(1/2) – 2*M*r + r^2))/(r^2*(M*r)^(1/2)*(2*a*(M*r)^(1/2) – 3*m*r + r^2)^2)

I tried: [Sr] = solve(D2,r) and got after some long warning that is pointless to paste here:

Sr =

0.00001*z1^2

So, where did the parameter z1 come from, and how is it defined? Was it even able to solve the problem?

Thx in advance!

Best Answer

You are not going to get anything useful.

Sr = RootOf(24*z^18 + (4*m*M^2 - 216*m*M - 16*M^2)*z^16 + 144*z^15*M^2*a + (-36*M^3*m^2 + 648*M^2*m^2 + 144*m*M^3)*z^14 + (27*M^4*m*a - 96*M^4*a - 864*M^3*a*m - 2*M^3*a*m^2)*z^13 + (-648*m^3*M^3 - 432*m^2*M^4 + 288*M^4*a^2 + 108*M^4*m^3)*z^12 + (-219*m^2*a*M^5 + 576*m*a*M^5 + 1296*m^2*M^4*a + 36*M^4*a*m^3 + 18*M^6*a*m)*z^11 + (-108*M^4*m^3 + 144*m^2*M^5 + 432*M^5*m^3 - 72*m^3*M^6 - 192*a^2*M^6 - 36*m^3*a^2*M^4 + 130*M^6*a^2*m - 12*M^5*a^2*m^2 - 864*a^2*m*M^5 - 48*m*M^6)*z^10 + (-90*M^6*a*m + 135*m^2*a*M^5 + 75*M^5*a^3*m^2 + 198*M^7*a*m^2 - 864*m^2*a*M^6 - 90*M^6*a^3*m - 126*M^6*a*m^3 + 192*M^6*a^3 + 54*M^5*a*m^4)*z^9 + (420*M^7*a^2*m^2 + 576*M^7*a^2*m - 36*M^6*a^2*m - 372*M^8*a^2*m - 108*m^3*M^6*a^2)*z^8 + (-216*M^7*a*m^2 + 27*M^6*a*m^3 - 51*M^6*a^3*m^3 + 132*M^8*a*m + 172*M^8*a^3*m + 36*m^4*a*M^7 - 124*m^2*a^3*M^7 - 128*M^8*a^3 - 36*M^8*a*m^3)*z^7 + (54*M^7*a^2*m^2 + 648*M^9*a^2*m^2 + 36*M^8*a^2*m - 60*M^8*a^4*m + 66*M^7*a^4*m^2 - 648*M^8*a^2*m^3)*z^6 + (-144*M^10*m*a^3 - 27*M^8*a^3*m + 918*m^3*M^8*a^3 + 54*M^7*a^3*m^4 + 162*M^8*a*m^3 - 108*M^9*a*m^2 - 828*M^9*a^3*m^2)*z^5 + (4*M^9*a^4*m^2 + 24*M^10*a^2*m - 108*M^8*a^2*m^3 - 36*M^9*a^2*m^2 - 396*M^8*a^4*m^3 - 576*M^10*m^3*a^2 + 576*M^9*a^2*m^4 + 392*M^10*a^4*m)*z^4 + (-1080*M^9*a^3*m^4 + 180*M^9*m^2*a^5 + 432*m^2*M^11*a^3 - 180*M^10*m*a^5 + 81*M^9*a^3*m^2 + 648*M^10*m^3*a^3)*z^3 + (-18*M^10*a^4*m - 912*m^2*M^11*a^4 + 576*M^10*a^4*m^3 - 96*m*M^12*a^4 + 432*M^9*a^4*m^4)*z^2 + (256*m*M^12*a^5 + 188*M^11*m^2*a^5 - 444*M^10*a^5*m^3)*z + 120*M^11*a^6*m^2 - 120*m*M^12*a^6, z)^2/M

The notation RootOf(expression, z) means the set of z such that the expression evalutes to 0 -- that is, the roots of the equation.

In this case you have a polynomial of order 18, and such polynomials (usually) do not have closed algebraic solutions.

Related Solutions

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.

MATLAB: Is this variable “undefined”

Your loop

for n1 = [2:.01:1.5];

has its end point (1.5) greater than its start point (2), so it is never executed.

If you want a loop to count down, you must use a negative increment (e.g., -0.01) rather than a positive one.

Best Answer

Related Solutions

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

MATLAB: Is this variable “undefined”

Related Question