problem with solve y sym

Question

javiera de la carrera 2017-10-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/362335-problem-with-solve-y-sym

评论： Walter Roberson 2017-10-23

Hi, I am trying to solve a system of 11 equations with symbolic variables:

S=[q1 - q2.*(0.11.*x1aux) - q4.*(sumxxaux+0.17.*x2aux) - fc1aux + vpasaj1aux - lf1 - v1.*Delayxaux - q2*0.11.*x1aux - x1aux.*v1.*exp(0.3487*sumq0aux/6).*0.3487./6==0;
   q6 - q7.*(0.17.*x2aux) - q4.*(sumxxaux+0.11.*x1aux) - fc2aux + vpasaj2aux - lf2 - v2.*Delayxaux - q7*0.17.*x2aux - x2aux.*v2.*exp(0.3487*sumq0aux/6).*0.3487./6==0;
   q4/q2-dif2==0];
P= [q51 - q31.*xaux(1,1) - q4.*(0.11.*x1aux+0.17.*x2aux+sumxxaux-xaux(1,1)) - fcaux(1,1) + vpasajaux(1,1) - lfaux(1,1) - vaux(1,1).*Delayxaux - q31.*xaux(1,1) - xaux(1,1).*vaux(1,1).*exp(0.3487*sumq0aux/6).*0.3487./6==0;
    q52 - q32.*xaux(1,2) - q4.*(0.11.*x1aux+0.17.*x2aux+sumxxaux-xaux(1,2)) - fcaux(1,2) + vpasajaux(1,2) - lfaux(1,2) - vaux(1,2).*Delayxaux - q32.*xaux(1,2) - xaux(1,2).*vaux(1,2).*exp(0.3487*sumq0aux/6).*0.3487./6==0;
    q53 - q33.*xaux(1,3) - q4.*(0.11.*x1aux+0.17.*x2aux+sumxxaux-xaux(1,3)) - fcaux(1,3) + vpasajaux(1,3) - lfaux(1,3) - vaux(1,3).*Delayxaux - q33.*xaux(1,3) - xaux(1,3).*vaux(1,3).*exp(0.3487*sumq0aux/6).*0.3487./6==0];
R= [(-(3*q4^4 - 2*q4^3*q31 - 2*q4^3*q32 - 2*q4^3*q33 - 2*q4^3*q7 + q4^2*q7*q31 + q4^2*q7*q32 + q4^2*q7*q33 + q4^2*q31*q32 + q4^2*q31*q33 + q4^2*q32*q33 - q7*q31*q32*q33)/(3*q2*q4^4 + 3*q4^4*q7 + 3*q4^4*q31 + 3*q4^4*q32 + 3*q4^4*q33 - 4*q4^5 - 2*q2*q4^3*q7 - 2*q2*q4^3*q31 - 2*q2*q4^3*q32 - 2*q2*q4^3*q33 - 2*q4^3*q7*q31 - 2*q4^3*q7*q32 - 2*q4^3*q7*q33 - 2*q4^3*q31*q32 - 2*q4^3*q31*q33 - 2*q4^3*q32*q33 + q2*q4^2*q7*q31 + q2*q4^2*q7*q32 + q2*q4^2*q7*q33 + q2*q4^2*q31*q32 + q2*q4^2*q31*q33 + q2*q4^2*q32*q33 + q4^2*q7*q31*q32 + q4^2*q7*q31*q33 + q4^2*q7*q32*q33 + q4^2*q31*q32*q33 - q2*q7*q31*q32*q33))-(eta.*(0.11.*x1aux))./(q1 - q2.*(0.11.*x1aux) - q4.*(sumxxaux+0.17.*x2aux))==0;
    (-(3*q4^4 - 2*q4^3*q31 - 2*q4^3*q32 - 2*q4^3*q33 - 2*q2*q4^3 + q2*q4^2*q31 + q2*q4^2*q32 + q2*q4^2*q33 + q4^2*q31*q32 + q4^2*q31*q33 + q4^2*q32*q33 - q2*q31*q32*q33)/(3*q2*q4^4 + 3*q4^4*q7 + 3*q4^4*q31 + 3*q4^4*q32 + 3*q4^4*q33 - 4*q4^5 - 2*q2*q4^3*q7 - 2*q2*q4^3*q31 - 2*q2*q4^3*q32 - 2*q2*q4^3*q33 - 2*q4^3*q7*q31 - 2*q4^3*q7*q32 - 2*q4^3*q7*q33 - 2*q4^3*q31*q32 - 2*q4^3*q31*q33 - 2*q4^3*q32*q33 + q2*q4^2*q7*q31 + q2*q4^2*q7*q32 + q2*q4^2*q7*q33 + q2*q4^2*q31*q32 + q2*q4^2*q31*q33 + q2*q4^2*q32*q33 + q4^2*q7*q31*q32 + q4^2*q7*q31*q33 + q4^2*q7*q32*q33 + q4^2*q31*q32*q33 - q2*q7*q31*q32*q33))-(eta.*(0.17.*x2aux))./(q6 - q7.*(0.17.*x2aux) - q4.*(0.11.*x1aux+sumxxaux))==0];
G= [(-(3*q4^4 - 2*q4^3*q7 - 2*q4^3*q32 - 2*q4^3*q33 - 2*q2*q4^3 + q2*q4^2*q7 + q2*q4^2*q32 + q2*q4^2*q33 + q4^2*q7*q32 + q4^2*q7*q33 + q4^2*q32*q33 - q2*q7*q32*q33)/(3*q2*q4^4 + 3*q4^4*q7 + 3*q4^4*q31 + 3*q4^4*q32 + 3*q4^4*q33 - 4*q4^5 - 2*q2*q4^3*q7 - 2*q2*q4^3*q31 - 2*q2*q4^3*q32 - 2*q2*q4^3*q33 - 2*q4^3*q7*q31 - 2*q4^3*q7*q32 - 2*q4^3*q7*q33 - 2*q4^3*q31*q32 - 2*q4^3*q31*q33 - 2*q4^3*q32*q33 + q2*q4^2*q7*q31 + q2*q4^2*q7*q32 + q2*q4^2*q7*q33 + q2*q4^2*q31*q32 + q2*q4^2*q31*q33 + q2*q4^2*q32*q33 + q4^2*q7*q31*q32 + q4^2*q7*q31*q33 + q4^2*q7*q32*q33 + q4^2*q31*q32*q33 - q2*q7*q31*q32*q33))-(eta.*(xaux(1,1)))./(q51 - q31.*xaux(1,1) - q4.*(0.11.*x1aux+0.17.*x2aux+sumxxaux-xaux(1,1)))==0;
    (-(3*q4^4 - 2*q4^3*q7 - 2*q4^3*q31 - 2*q4^3*q33 - 2*q2*q4^3 + q2*q4^2*q7 + q2*q4^2*q31 + q2*q4^2*q33 + q4^2*q7*q31 + q4^2*q7*q33 + q4^2*q31*q33 - q2*q7*q31*q33)/(3*q2*q4^4 + 3*q4^4*q7 + 3*q4^4*q31 + 3*q4^4*q32 + 3*q4^4*q33 - 4*q4^5 - 2*q2*q4^3*q7 - 2*q2*q4^3*q31 - 2*q2*q4^3*q32 - 2*q2*q4^3*q33 - 2*q4^3*q7*q31 - 2*q4^3*q7*q32 - 2*q4^3*q7*q33 - 2*q4^3*q31*q32 - 2*q4^3*q31*q33 - 2*q4^3*q32*q33 + q2*q4^2*q7*q31 + q2*q4^2*q7*q32 + q2*q4^2*q7*q33 + q2*q4^2*q31*q32 + q2*q4^2*q31*q33 + q2*q4^2*q32*q33 + q4^2*q7*q31*q32 + q4^2*q7*q31*q33 + q4^2*q7*q32*q33 + q4^2*q31*q32*q33 - q2*q7*q31*q32*q33))-(eta.*(xaux(1,2)))./(q52 - q32.*xaux(1,2) - q4.*(0.11.*x1aux+0.17.*x2aux+sumxxaux-xaux(1,2)))==0;
    (-(3*q4^4 - 2*q4^3*q7 - 2*q4^3*q31 - 2*q4^3*q32 - 2*q2*q4^3 + q2*q4^2*q7 + q2*q4^2*q31 + q2*q4^2*q32 + q4^2*q7*q31 + q4^2*q7*q32 + q4^2*q31*q32 - q2*q7*q31*q32)/(3*q2*q4^4 + 3*q4^4*q7 + 3*q4^4*q31 + 3*q4^4*q32 + 3*q4^4*q33 - 4*q4^5 - 2*q2*q4^3*q7 - 2*q2*q4^3*q31 - 2*q2*q4^3*q32 - 2*q2*q4^3*q33 - 2*q4^3*q7*q31 - 2*q4^3*q7*q32 - 2*q4^3*q7*q33 - 2*q4^3*q31*q32 - 2*q4^3*q31*q33 - 2*q4^3*q32*q33 + q2*q4^2*q7*q31 + q2*q4^2*q7*q32 + q2*q4^2*q7*q33 + q2*q4^2*q31*q32 + q2*q4^2*q31*q33 + q2*q4^2*q32*q33 + q4^2*q7*q31*q32 + q4^2*q7*q31*q33 + q4^2*q7*q32*q33 + q4^2*q31*q32*q33 - q2*q7*q31*q32*q33))-(eta.*(xaux(1,3)))./(q53 - q33.*xaux(1,3) - q4.*(0.11.*x1aux+0.17.*x2aux+sumxxaux-xaux(1,3)))==0];
D=[S;R;P;G];
[t11,t22,t331,t332,t333,t44,t551,t552,t553,t66,t77,param,cond]=solve(D,q1,q2,q31,q32,q33,q4,q51,q52,q53,q6,q7,'ReturnConditions',true,'Real',true);

But I obtain this answer:

Error using mupadengine/feval (line 163)
The row index is out of range.
Error in solve (line 369)
            varargout{i} = transpose(eng.feval('map', solutions, '_index', i));

I don't know why! please help me!

Thanks

6 个评论
显示 4更早的评论隐藏 4更早的评论

Walter Roberson 2017-10-20

编辑：Walter Roberson 2017-10-20

在 MATLAB Online 中打开

It looks to me as if the code could be tested if it were preceded with

[Delayxaux, eta, fc1aux, fc2aux, lf1, lf2, sumxxaux, sumq0aux, v1, v2, vpasaj1aux, vpasaj2aux, x1aux, x2aux] = deal(rand, rand, rand, rand, rand, rand, rand, rand, rand, rand, rand, rand, rand, rand);
fcaux = rand(1,3);
lfaux = rand(1,3);
vaux = rand(1,3);
vpasajaux = rand(1,3);
xaux = rand(1,3);
syms dif2 q1 q2 q4 q6 q7 q31 q32 q33 q51 q52 q53

You appear to be solving 11 equations for 11 out of the 12 variables, so you would expect most of the outputs to be in terms of dif2

Can any constraints be put on the symbols? Such as real valued, or non-negative?

javiera de la carrera 2017-10-20

Yes, this is the code: syms q1 positive

syms q2 positive

syms q31 positive

syms q32 positive

syms q33 positive

syms q4 positive

syms q51 positive

syms q52 positive

syms q53 positive

syms q6 positive %a united

syms q7 positive %b united

syms dif2

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2017-10-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/362335-problem-with-solve-y-sym#answer_286885

在 MATLAB Online 中打开

With the above initialization, in R2017b, it works hard on it for a while and then says "Warning: Unable to find explicit solution." and returns empty symbols.

This does not surprise me. If you

collect(D, q4)

then you will see that you have q4^5 in some of the equations. There is no general closed form solution to equations of degree 5 and it is not possible for MATLAB to look for numeric approximations because of the free variable dif2 .

5 个评论
显示 3更早的评论隐藏 3更早的评论

Walter Roberson 2017-10-21

编辑：Walter Roberson 2017-10-21

在 MATLAB Online 中打开

Cleaner version

[Delayxaux, eta, fc1aux, fc2aux, lf1, lf2, sumxxaux, sumq0aux, v1, v2, vpasaj1aux, vpasaj2aux, x1aux, x2aux] = deal(rand, rand, rand, rand, rand, rand, rand, rand, rand, rand, rand, rand, rand, rand);
fcaux = rand(1,3);
lfaux = rand(1,3);
vaux = rand(1,3);
vpasajaux = rand(1,3);
xaux = rand(1,3);
syms dif2 
syms q1 q2 q4 q6 q7 q31 q32 q33 q51 q52 q53 positive
S=[q1 - q2.*(0.11.*x1aux) - q4.*(sumxxaux+0.17.*x2aux) - fc1aux + vpasaj1aux - lf1 - v1.*Delayxaux - q2*0.11.*x1aux - x1aux.*v1.*exp(0.3487*sumq0aux/6).*0.3487./6==0;
   q6 - q7.*(0.17.*x2aux) - q4.*(sumxxaux+0.11.*x1aux) - fc2aux + vpasaj2aux - lf2 - v2.*Delayxaux - q7*0.17.*x2aux - x2aux.*v2.*exp(0.3487*sumq0aux/6).*0.3487./6==0;
   q4/q2-dif2==0];
P= [q51 - q31.*xaux(1,1) - q4.*(0.11.*x1aux+0.17.*x2aux+sumxxaux-xaux(1,1)) - fcaux(1,1) + vpasajaux(1,1) - lfaux(1,1) - vaux(1,1).*Delayxaux - q31.*xaux(1,1) - xaux(1,1).*vaux(1,1).*exp(0.3487*sumq0aux/6).*0.3487./6==0;
    q52 - q32.*xaux(1,2) - q4.*(0.11.*x1aux+0.17.*x2aux+sumxxaux-xaux(1,2)) - fcaux(1,2) + vpasajaux(1,2) - lfaux(1,2) - vaux(1,2).*Delayxaux - q32.*xaux(1,2) - xaux(1,2).*vaux(1,2).*exp(0.3487*sumq0aux/6).*0.3487./6==0;
    q53 - q33.*xaux(1,3) - q4.*(0.11.*x1aux+0.17.*x2aux+sumxxaux-xaux(1,3)) - fcaux(1,3) + vpasajaux(1,3) - lfaux(1,3) - vaux(1,3).*Delayxaux - q33.*xaux(1,3) - xaux(1,3).*vaux(1,3).*exp(0.3487*sumq0aux/6).*0.3487./6==0];
R= [(-(3*q4^4 - 2*q4^3*q31 - 2*q4^3*q32 - 2*q4^3*q33 - 2*q4^3*q7 + q4^2*q7*q31 + q4^2*q7*q32 + q4^2*q7*q33 + q4^2*q31*q32 + q4^2*q31*q33 + q4^2*q32*q33 - q7*q31*q32*q33)/(3*q2*q4^4 + 3*q4^4*q7 + 3*q4^4*q31 + 3*q4^4*q32 + 3*q4^4*q33 - 4*q4^5 - 2*q2*q4^3*q7 - 2*q2*q4^3*q31 - 2*q2*q4^3*q32 - 2*q2*q4^3*q33 - 2*q4^3*q7*q31 - 2*q4^3*q7*q32 - 2*q4^3*q7*q33 - 2*q4^3*q31*q32 - 2*q4^3*q31*q33 - 2*q4^3*q32*q33 + q2*q4^2*q7*q31 + q2*q4^2*q7*q32 + q2*q4^2*q7*q33 + q2*q4^2*q31*q32 + q2*q4^2*q31*q33 + q2*q4^2*q32*q33 + q4^2*q7*q31*q32 + q4^2*q7*q31*q33 + q4^2*q7*q32*q33 + q4^2*q31*q32*q33 - q2*q7*q31*q32*q33))-(eta.*(0.11.*x1aux))./(q1 - q2.*(0.11.*x1aux) - q4.*(sumxxaux+0.17.*x2aux))==0;
    (-(3*q4^4 - 2*q4^3*q31 - 2*q4^3*q32 - 2*q4^3*q33 - 2*q2*q4^3 + q2*q4^2*q31 + q2*q4^2*q32 + q2*q4^2*q33 + q4^2*q31*q32 + q4^2*q31*q33 + q4^2*q32*q33 - q2*q31*q32*q33)/(3*q2*q4^4 + 3*q4^4*q7 + 3*q4^4*q31 + 3*q4^4*q32 + 3*q4^4*q33 - 4*q4^5 - 2*q2*q4^3*q7 - 2*q2*q4^3*q31 - 2*q2*q4^3*q32 - 2*q2*q4^3*q33 - 2*q4^3*q7*q31 - 2*q4^3*q7*q32 - 2*q4^3*q7*q33 - 2*q4^3*q31*q32 - 2*q4^3*q31*q33 - 2*q4^3*q32*q33 + q2*q4^2*q7*q31 + q2*q4^2*q7*q32 + q2*q4^2*q7*q33 + q2*q4^2*q31*q32 + q2*q4^2*q31*q33 + q2*q4^2*q32*q33 + q4^2*q7*q31*q32 + q4^2*q7*q31*q33 + q4^2*q7*q32*q33 + q4^2*q31*q32*q33 - q2*q7*q31*q32*q33))-(eta.*(0.17.*x2aux))./(q6 - q7.*(0.17.*x2aux) - q4.*(0.11.*x1aux+sumxxaux))==0];
G= [(-(3*q4^4 - 2*q4^3*q7 - 2*q4^3*q32 - 2*q4^3*q33 - 2*q2*q4^3 + q2*q4^2*q7 + q2*q4^2*q32 + q2*q4^2*q33 + q4^2*q7*q32 + q4^2*q7*q33 + q4^2*q32*q33 - q2*q7*q32*q33)/(3*q2*q4^4 + 3*q4^4*q7 + 3*q4^4*q31 + 3*q4^4*q32 + 3*q4^4*q33 - 4*q4^5 - 2*q2*q4^3*q7 - 2*q2*q4^3*q31 - 2*q2*q4^3*q32 - 2*q2*q4^3*q33 - 2*q4^3*q7*q31 - 2*q4^3*q7*q32 - 2*q4^3*q7*q33 - 2*q4^3*q31*q32 - 2*q4^3*q31*q33 - 2*q4^3*q32*q33 + q2*q4^2*q7*q31 + q2*q4^2*q7*q32 + q2*q4^2*q7*q33 + q2*q4^2*q31*q32 + q2*q4^2*q31*q33 + q2*q4^2*q32*q33 + q4^2*q7*q31*q32 + q4^2*q7*q31*q33 + q4^2*q7*q32*q33 + q4^2*q31*q32*q33 - q2*q7*q31*q32*q33))-(eta.*(xaux(1,1)))./(q51 - q31.*xaux(1,1) - q4.*(0.11.*x1aux+0.17.*x2aux+sumxxaux-xaux(1,1)))==0;
    (-(3*q4^4 - 2*q4^3*q7 - 2*q4^3*q31 - 2*q4^3*q33 - 2*q2*q4^3 + q2*q4^2*q7 + q2*q4^2*q31 + q2*q4^2*q33 + q4^2*q7*q31 + q4^2*q7*q33 + q4^2*q31*q33 - q2*q7*q31*q33)/(3*q2*q4^4 + 3*q4^4*q7 + 3*q4^4*q31 + 3*q4^4*q32 + 3*q4^4*q33 - 4*q4^5 - 2*q2*q4^3*q7 - 2*q2*q4^3*q31 - 2*q2*q4^3*q32 - 2*q2*q4^3*q33 - 2*q4^3*q7*q31 - 2*q4^3*q7*q32 - 2*q4^3*q7*q33 - 2*q4^3*q31*q32 - 2*q4^3*q31*q33 - 2*q4^3*q32*q33 + q2*q4^2*q7*q31 + q2*q4^2*q7*q32 + q2*q4^2*q7*q33 + q2*q4^2*q31*q32 + q2*q4^2*q31*q33 + q2*q4^2*q32*q33 + q4^2*q7*q31*q32 + q4^2*q7*q31*q33 + q4^2*q7*q32*q33 + q4^2*q31*q32*q33 - q2*q7*q31*q32*q33))-(eta.*(xaux(1,2)))./(q52 - q32.*xaux(1,2) - q4.*(0.11.*x1aux+0.17.*x2aux+sumxxaux-xaux(1,2)))==0;
    (-(3*q4^4 - 2*q4^3*q7 - 2*q4^3*q31 - 2*q4^3*q32 - 2*q2*q4^3 + q2*q4^2*q7 + q2*q4^2*q31 + q2*q4^2*q32 + q4^2*q7*q31 + q4^2*q7*q32 + q4^2*q31*q32 - q2*q7*q31*q32)/(3*q2*q4^4 + 3*q4^4*q7 + 3*q4^4*q31 + 3*q4^4*q32 + 3*q4^4*q33 - 4*q4^5 - 2*q2*q4^3*q7 - 2*q2*q4^3*q31 - 2*q2*q4^3*q32 - 2*q2*q4^3*q33 - 2*q4^3*q7*q31 - 2*q4^3*q7*q32 - 2*q4^3*q7*q33 - 2*q4^3*q31*q32 - 2*q4^3*q31*q33 - 2*q4^3*q32*q33 + q2*q4^2*q7*q31 + q2*q4^2*q7*q32 + q2*q4^2*q7*q33 + q2*q4^2*q31*q32 + q2*q4^2*q31*q33 + q2*q4^2*q32*q33 + q4^2*q7*q31*q32 + q4^2*q7*q31*q33 + q4^2*q7*q32*q33 + q4^2*q31*q32*q33 - q2*q7*q31*q32*q33))-(eta.*(xaux(1,3)))./(q53 - q33.*xaux(1,3) - q4.*(0.11.*x1aux+0.17.*x2aux+sumxxaux-xaux(1,3)))==0];
D=[S;R;P;G];
%[t11,t22,t331,t332,t333,t44,t551,t552,t553,t66,t77,param,cond]=solve(D,q1,q2,q31,q32,q33,q4,q51,q52,q53,q6,q7,'ReturnConditions',true,'Real',true);
%sol = solve(D,q1,q2,q31,q32,q33,q4,q51,q52,q53,q6,q7,'ReturnConditions',true);
Q4 = solve(D(1), q4, 'ReturnConditions', true)
D2 = subs(D, q4, Q4.q4)
Q1 = solve(D2(2), q1, 'ReturnConditions', true)
D3 = subs(D2, q1, Q1.q1)
Q2 = solve(D3(3), q2, 'ReturnConditions', true)
D4 = subs(D3, q2, Q2.q2)
Q31 = solve(D4(6), q31, 'ReturnConditions', true)
D5 = subs(D4, q31, Q31.q31)
Q32 = solve(D5(7), q32, 'ReturnConditions', true)
D6 = subs(D5, q32, Q32.q32)
Q33 = solve(D6(8), q33, 'ReturnConditions', true)
D7 = subs(D6, q33, Q33.q33)
Q51 = solve(D7(9), q51, 'ReturnConditions', true)
D8 = subs(D7, q51, Q51.q51)
%this next takes a while but does finish
Q52 = solve(D8(11), q52, 'ReturnConditions', true)
keyboard   %do not accidentally continue to the next one!
%gets really really bogged down in the following
D9 = subs(D8, q52, Q52.q52(1))

javiera de la carrera 2017-10-23

Thanks! I think I have to change the problem. Or, Do you think another solution? Like a command to approximate equations or simplify to the solve comand? I tried with simplify D but nothing change.

Thanks!

Walter Roberson 2017-10-23

在 MATLAB Online 中打开

Your R and your G both contain the ratio of quartics (polynomials of degree 4.) Just knowing that, we can predict that the solutions are going to be at least degree 7, or degree 8, or degree 16. You have 5 like that, and MATLAB needs to explore all of the possibilities for all of them. That is going to be slow and memory intensive.

Perhaps you can taylor with respect to q4 -- but to what order?

I extracted the division subexpression from the first R entry, and tossed in some random values for some of the variables

t2 = (3*q4^4 - 2*q4^3*q31 - 2*q4^3*q32 - 2*q4^3*q33 - 2*q4^3*q7 + q4^2*q7*q31 + q4^2*q7*q32 + q4^2*q7*q33 + q4^2*q31*q32 + q4^2*q31*q33 + q4^2*q32*q33 - q7*q31*q32*q33)/(3*q2*q4^4 + 3*q4^4*q7 + 3*q4^4*q31 + 3*q4^4*q32 + 3*q4^4*q33 - 4*q4^5 - 2*q2*q4^3*q7 - 2*q2*q4^3*q31 - 2*q2*q4^3*q32 - 2*q2*q4^3*q33 - 2*q4^3*q7*q31 - 2*q4^3*q7*q32 - 2*q4^3*q7*q33 - 2*q4^3*q31*q32 - 2*q4^3*q31*q33 - 2*q4^3*q32*q33 + q2*q4^2*q7*q31 + q2*q4^2*q7*q32 + q2*q4^2*q7*q33 + q2*q4^2*q31*q32 + q2*q4^2*q31*q33 + q2*q4^2*q32*q33 + q4^2*q7*q31*q32 + q4^2*q7*q31*q33 + q4^2*q7*q32*q33 + q4^2*q31*q32*q33 - q2*q7*q31*q32*q33)
t3 = subs(t2,[q2,q32,q33,q7],[2,32,33,7])
fsurf(t3,[-100 100 -80 40])

The result has at least 4 curved lines of discontinuity, indicating that if you cannot promise that the variables will stay in the "safe" ranges, that there cannot be any finite polynomial approximation.

I don't know... perhaps you could do something with continued fractions.

请先登录，再进行评论。