Error using symengine division by zero

I am trying to solve a vibration by using the Rayleigh-Ritz method with an admissable function. The error was as follows - Error using symengine Division by zero. Error in sym/subs>mupadsubs ...

WONG Weng Tuck

I am trying to solve a vibration by using the Rayleigh-Ritz method with an admissable function.

The error was as follows —

Error using symengine

Division by zero.

Error in sym/subs>mupadsubs (line 160)

G = mupadmex(‘symobj::fullsubs’,F.s,X2,Y2);

Error in sym/subs (line 145)

G = mupadsubs(F,X,Y);

Error in ca2draft3 (line 54)

MassMatrix(rows, columns) = subs(Mij, [i,j], [rows, columns])

Tried using eps and still got the error. Could anyone please assist?

Code is as follows —

AreaWing = 0.3 — (y/440) + ((y^2)/9680);

SecMomAreaWing = 0.001 — (y/132000) + ((y^2)/3872000) — ((y^3)/85184000);

phi_i_y = ((y^2)/(44^2))*cos((i*pi*y)/44);

phi_j_y = ((y^2)/(44^2))*cos((j*pi*y)/44);

phi_doubleprime_i_y = diff(phi_i_y,y,2);

phi_doubleprime_j_y = diff(phi_j_y,y,2);

phi_i_inbdEng = ((8.8^2)/44^2)*cos(0.2*i*pi);

phi_i_midEng = ((17.6^2)/44^2)*cos(0.4*i*pi);

phi_i_outerEng = ((26.4^2)/44^2)*cos(0.6*i*pi);

phi_j_inbdEng = ((8.8^2)/44^2)*cos(0.2*j*pi);

phi_j_midEng = ((17.6^2)/44^2)*cos(0.4*j*pi);

phi_j_outerEng = ((26.4^2)/44^2)*cos(0.6*j*pi);

MassinbdEng = 4100*(phi_i_inbdEng)*(phi_j_inbdEng);

MassmidEng = 4100*(phi_i_midEng)*(phi_j_midEng);

MassouterEng = 4100*(phi_i_outerEng)*(phi_j_outerEng);

StiffnessWing = 2e12*(SecMomAreaWing)*(phi_doubleprime_i_y)*(phi_doubleprime_j_y);

MassWing = 300*(AreaWing)*(phi_i_y)*(phi_j_y);

Mij = int(MassWing,y,0,44) + MassinbdEng + MassmidEng + MassouterEng;

Kij = int(StiffnessWing,y,0,44);

MassMatrix(rows, columns) = subs(Mij, [i,j], [rows, columns])

StiffMatrix(rows, columns) = subs(Kij, [i,j], [rows, columns])

Systemmatrix = inv(MassMatrix)*StiffMatrix;

[V,D] = eig(Systemmatrix);

[D_sorted, ind] = sort(diag(D),‘ascend’);

nat_freq_one = sqrt(D_sorted(1))

nat_freq_two = sqrt(D_sorted(2))

mode_shape_one = V_sorted(:,1)

mode_shape_two = V_sorted(:,2)

nat_freq_three = sqrt(D_sorted(3))

nat_freq_four = sqrt(D_sorted(4))

mode_shape_three = V_sorted(:,3)

mode_shape_four = V_sorted(:,4)

Accepted Answer

Walter Roberson

That integral ends up with a division by (i-j) . When rows == cols that is a division by 0.

You can get around that by subs() for i, and then taking the limit to j

i, j, y are indefinite variables at that point. They do not «equal» anything at that point.

Furthermore, you later substitute in row and column numbers for i and j, so i and j are always positive.

Here is the new loop:

MassMatrix(rows, columns) = simplify(limit(subs(Mij, i, rows), j, columns));

StiffMatrix(rows, columns) = simplify(limit(subs(Kij, i, rows), j, columns));


More Answers (0)

See Also

Categories

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

An Error Occurred

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Источник

Fireman

Пользователь
Сообщения: 5
Зарегистрирован: Ср фев 06, 2019 1:12 am

упрощение выражений

Приветствую

Подскажите со следующим вопросом:
у меня есть определённое «многоэтажное» выражение, которое очень бы хотелось упростить
Изображение

Очень хотелось бы формулу вида a/b, т.е. двухуровневую (исходник 16уровневый) и очевидно что такая формула существует, пусть и длинная.

Если я использую коману simplify, то получаю 8уровневую упрощённу формулу, что меня не устраивает.

Подскажите как заставить mathlab сделать упрощённую 2уровневую формулу.

P.S.

сам формулу

Код: Выделить всё

t =

((((z_1 - z_2)/(x_1^2 - x_2^2) - (z_1 - z_3)/(x_1^2 - x_3^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_3)/(x_1^2 - x_3^2)) - ((z_1 - z_2)/(x_1^2 - x_2^2) - (z_1 - z_4)/(x_1^2 - x_4^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_4)/(x_1^2 - x_4^2)))/(((y_1^2 - y_2^2)/(x_1^2 - x_2^2) - (y_1^2 - y_3^2)/(x_1^2 - x_3^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_3)/(x_1^2 - x_3^2)) - ((y_1^2 - y_2^2)/(x_1^2 - x_2^2) - (y_1^2 - y_4^2)/(x_1^2 - x_4^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_4)/(x_1^2 - x_4^2))) - (((z_1 - z_2)/(x_1^2 - x_2^2) - (z_1 - z_3)/(x_1^2 - x_3^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_3)/(x_1^2 - x_3^2)) - ((z_1 - z_2)/(x_1^2 - x_2^2) - (z_1 - z_5)/(x_1^2 - x_5^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_5)/(x_1^2 - x_5^2)))/(((y_1^2 - y_2^2)/(x_1^2 - x_2^2) - (y_1^2 - y_3^2)/(x_1^2 - x_3^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_3)/(x_1^2 - x_3^2)) - ((y_1^2 - y_2^2)/(x_1^2 - x_2^2) - (y_1^2 - y_5^2)/(x_1^2 - x_5^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_5)/(x_1^2 - x_5^2))))/((((y_1 - y_2)/(x_1^2 - x_2^2) - (y_1 - y_3)/(x_1^2 - x_3^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_3)/(x_1^2 - x_3^2)) - ((y_1 - y_2)/(x_1^2 - x_2^2) - (y_1 - y_4)/(x_1^2 - x_4^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_4)/(x_1^2 - x_4^2)))/(((y_1^2 - y_2^2)/(x_1^2 - x_2^2) - (y_1^2 - y_3^2)/(x_1^2 - x_3^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_3)/(x_1^2 - x_3^2)) - ((y_1^2 - y_2^2)/(x_1^2 - x_2^2) - (y_1^2 - y_4^2)/(x_1^2 - x_4^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_4)/(x_1^2 - x_4^2))) - (((y_1 - y_2)/(x_1^2 - x_2^2) - (y_1 - y_3)/(x_1^2 - x_3^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_3)/(x_1^2 - x_3^2)) - ((y_1 - y_2)/(x_1^2 - x_2^2) - (y_1 - y_5)/(x_1^2 - x_5^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_5)/(x_1^2 - x_5^2)))/(((y_1^2 - y_2^2)/(x_1^2 - x_2^2) - (y_1^2 - y_3^2)/(x_1^2 - x_3^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_3)/(x_1^2 - x_3^2)) - ((y_1^2 - y_2^2)/(x_1^2 - x_2^2) - (y_1^2 - y_5^2)/(x_1^2 - x_5^2))/((x_1 - x_2)/(x_1^2 - x_2^2) - (x_1 - x_5)/(x_1^2 - x_5^2))))

выдается такой результат:

Код: Выделить всё

ans =

(((x_1^2 z_2 - x_2^2 z_1 - x_1^2 z_3 + x_3^2 z_1 + x_2^2 z_3 - x_3^2 z_2)/((x_1 - x_2) (x_1 - x_3) (x_2 - x_3)) - (x_1^2 z_2 - x_2^2 z_1 - x_1^2 z_4 + x_4^2 z_1 + x_2^2 z_4 - x_4^2 z_2)/((x_1 - x_2) (x_1 - x_4) (x_2 - x_4)))/((x_1^2 y_2^2 - x_1^2 y_3^2 - x_2^2 y_1^2 + x_2^2 y_3^2 + x_3^2 y_1^2 - x_3^2 y_2^2)/((x_1 - x_2) (x_1 - x_3) (x_2 - x_3)) - (x_1^2 y_2^2 - x_1^2 y_4^2 - x_2^2 y_1^2 + x_2^2 y_4^2 + x_4^2 y_1^2 - x_4^2 y_2^2)/((x_1 - x_2) (x_1 - x_4) (x_2 - x_4))) - ((x_1^2 z_2 - x_2^2 z_1 - x_1^2 z_3 + x_3^2 z_1 + x_2^2 z_3 - x_3^2 z_2)/((x_1 - x_2) (x_1 - x_3) (x_2 - x_3)) - (x_1^2 z_2 - x_2^2 z_1 - x_1^2 z_5 + x_5^2 z_1 + x_2^2 z_5 - x_5^2 z_2)/((x_1 - x_2) (x_1 - x_5) (x_2 - x_5)))/((x_1^2 y_2^2 - x_1^2 y_3^2 - x_2^2 y_1^2 + x_2^2 y_3^2 + x_3^2 y_1^2 - x_3^2 y_2^2)/((x_1 - x_2) (x_1 - x_3) (x_2 - x_3)) - (x_1^2 y_2^2 - x_1^2 y_5^2 - x_2^2 y_1^2 + x_2^2 y_5^2 + x_5^2 y_1^2 - x_5^2 y_2^2)/((x_1 - x_2) (x_1 - x_5) (x_2 - x_5))))/(((x_1^2 y_2 - x_2^2 y_1 - x_1^2 y_3 + x_3^2 y_1 + x_2^2 y_3 - x_3^2 y_2)/((x_1 - x_2) (x_1 - x_3) (x_2 - x_3)) - (x_1^2 y_2 - x_2^2 y_1 - x_1^2 y_4 + x_4^2 y_1 + x_2^2 y_4 - x_4^2 y_2)/((x_1 - x_2) (x_1 - x_4) (x_2 - x_4)))/((x_1^2 y_2^2 - x_1^2 y_3^2 - x_2^2 y_1^2 + x_2^2 y_3^2 + x_3^2 y_1^2 - x_3^2 y_2^2)/((x_1 - x_2) (x_1 - x_3) (x_2 - x_3)) - (x_1^2 y_2^2 - x_1^2 y_4^2 - x_2^2 y_1^2 + x_2^2 y_4^2 + x_4^2 y_1^2 - x_4^2 y_2^2)/((x_1 - x_2) (x_1 - x_4) (x_2 - x_4))) - ((x_1^2 y_2 - x_2^2 y_1 - x_1^2 y_3 + x_3^2 y_1 + x_2^2 y_3 - x_3^2 y_2)/((x_1 - x_2) (x_1 - x_3) (x_2 - x_3)) - (x_1^2 y_2 - x_2^2 y_1 - x_1^2 y_5 + x_5^2 y_1 + x_2^2 y_5 - x_5^2 y_2)/((x_1 - x_2) (x_1 - x_5) (x_2 - x_5)))/((x_1^2 y_2^2 - x_1^2 y_3^2 - x_2^2 y_1^2 + x_2^2 y_3^2 + x_3^2 y_1^2 - x_3^2 y_2^2)/((x_1 - x_2) (x_1 - x_3) (x_2 - x_3)) - (x_1^2 y_2^2 - x_1^2 y_5^2 - x_2^2 y_1^2 + x_2^2 y_5^2 + x_5^2 y_1^2 - x_5^2 y_2^2)/((x_1 - x_2) (x_1 - x_5) (x_2 - x_5))))

sandy

Эксперт
Сообщения: 5601
Зарегистрирован: Ср сен 22, 2004 4:49 pm

Re: упрощение выражений

Сообщение sandy » Чт фев 07, 2019 12:17 am

Попробуйте функцию simplifyFraction

С уважением

Александр Сергиенко

Fireman

Пользователь
Сообщения: 5
Зарегистрирован: Ср фев 06, 2019 1:12 am

Re: упрощение выражений

Сообщение Fireman » Чт фев 07, 2019 1:04 pm

Пока в отъезде, в субботу опробуют все предложенные способы

Подумалось — а не сработает ли многошаговое упрощение, типа

?

Fireman

Пользователь
Сообщения: 5
Зарегистрирован: Ср фев 06, 2019 1:12 am

Re: упрощение выражений

Сообщение Fireman » Сб фев 09, 2019 10:54 pm

sandy писал(а):Попробуйте функцию simplifyFraction

Огромное спасибо, сработало
Получилось выражение вида a/b — как раз то, что нужно

А не подскажешь со второй проблемой:
Я хочу произвести замены

Код: Выделить всё

syms d_x
syms d_y

x_2 = x_1 - d_x
y_2 = y_1

x_3 = x_1 + d_x
y_3 = y_1

x_4 = x_1
y_4 = y_1 - d_y

x_5 = x_1
y_5 = y_1 + d_y

в результате вместо результата выдается NaN

Если же я беру отдельно a и b из предыдущего результата и применяю к ним такие преобразования — то никаких Nan не возникает и получается корректное выражение

А потом я делаю a/b руками и опять корректное выражение вылезает
Откуда же тогда NaN изначально?

barkot

Пользователь
Сообщения: 61
Зарегистрирован: Вс янв 14, 2018 2:28 pm

Re: упрощение выражений

Сообщение barkot » Вс фев 10, 2019 9:16 am

а если так ?

Код: Выделить всё

t1=simplifyFraction(t);
subs(t1,{x2,y2,x3,y3,x4,y4,x5,y5},[x1-dx,y1,x1+dx,y1,x1,y1-dy,x1,y1+dy])

Fireman

Пользователь
Сообщения: 5
Зарегистрирован: Ср фев 06, 2019 1:12 am

Re: упрощение выражений

Сообщение Fireman » Пн фев 11, 2019 1:10 am

barkot писал(а):а если так ?

Код: Выделить всё

t1=simplifyFraction(t);
subs(t1,{x2,y2,x3,y3,x4,y4,x5,y5},[x1-dx,y1,x1+dx,y1,x1,y1-dy,x1,y1+dy])

Замечательно работает, но если делать так

Код: Выделить всё

subs(t,{x2,y2,x3,y3,x4,y4,x5,y5},[x1-dx,y1,x1+dx,y1,x1,y1-dy,x1,y1+dy])

т.е. если не делать начального упрощения, то выдает

Error using symengine
Division by zero.

Error in sym/subs>mupadsubs (line 160)
G = mupadmex(‘symobj::fullsubs’,F.s,X2,Y2);

Error in sym/subs (line 145)
G = mupadsubs(F,X,Y);

Error in script (line 43)
res = subs(t, {x_2, y_2, x_3, y_3, x_4, y_4, x_5, y_5}, [x_1 — d_x, y_1, x_1 + d_x, y_1, x_1, y_1 — d_y, x_1, y_1 + d_y]);

sandy

Эксперт
Сообщения: 5601
Зарегистрирован: Ср сен 22, 2004 4:49 pm

Re: упрощение выражений

Сообщение sandy » Пн фев 11, 2019 8:15 pm

Ну а чего вы хотите? Если сделать эти подстановки в вашем исходном выражении, действительно получается деление на ноль.

С уважением

Александр Сергиенко

Fireman

Пользователь
Сообщения: 5
Зарегистрирован: Ср фев 06, 2019 1:12 am

Re: упрощение выражений

Сообщение Fireman » Пн фев 11, 2019 11:06 pm

sandy писал(а):Ну а чего вы хотите? Если сделать эти подстановки в вашем исходном выражении, действительно получается деление на ноль.

кхм, действительно так, но это в общей записи, при сокращении этот эффект должен исчезнуть

вообще в посте приведена задача нахождения коэффициентов параболоида по 5 точкам, понятно что есть частные случаи, где будет деление на 0, но в общем случае его быть не должно

а когда я привожу нижние замены — это 5 точек, расположенные на кресте, при разных z1..z5 там вообще деления на 0 не должно возникать

sandy

Эксперт
Сообщения: 5601
Зарегистрирован: Ср сен 22, 2004 4:49 pm

Re: упрощение выражений

Сообщение sandy » Пн фев 11, 2019 11:43 pm

Fireman писал(а):при сокращении этот эффект должен исчезнуть

Ну так он и исчезает — ведь его нет, когда вы выполняете подстановку в упрощенном выражении t1. Не понимаю, что вам не нравится.

С уважением

Александр Сергиенко

Источник
clear all

syms s w 

G = 1/((s)*(s+1)*(s+2)); %transfer function
G_w = subs(G,s,j*w); 
W= [-100:0.01:100]; %[min_range:step size:max_range]
nyq = eval(subs(G_w,w,W)); 

x = real(nyq)
y = imag(nyq)

plot(x,y)

Кажется, я не могу запустить этот код, и он продолжает отображать ошибку в строке 100 ++, где у меня меньше 20 строк.

Error using symengine (line 59)
Division by zero.
Error in sym/subs>mupadsubs (line 139)
G = mupadmex('symobj::fullsubs',F.s,X2,Y2);
Error in sym/subs (line 124)
    G = mupadsubs(F,X,Y);
Error in nyquist2 (line 8)
nyq = eval(subs(G_w,w,W)); %replace W with w in equation G_w

Это показанные ошибки, может ли мне в этом помочь специалист?

2 ответа

Лучший ответ

Ошибка возникает из-за того, что вы вычисляете G_w с использованием массива W, и этот массив содержит значение 0, приводящее к делению на ноль и, следовательно, к ошибке.

Ошибка при использовании symengine (строка 59)
Деление на ноль.

Что вы можете сделать, чтобы обойти это, — заменить 0 в W на eps.

% Replace zeros with epsilon
W(W == 0) = eps;

nyq = eval(subs(G_w,w,W)); 

x = real(nyq)
y = imag(nyq)

plot(x,y)

В качестве примечания: ошибка не связана с проблемой со строкой 100+ вашего кода, а скорее, трассировка стека указывает, что ошибка на самом деле возникла из внутри функция, которую вы вызываете

Трассировка стека упорядочена от места возникновения ошибки до кода, который вы вызвали для ее создания.

Error using symengine (line 59)             <--- WHERE THE ERROR HAPPENED
Division by zero.                           <--- THIS IS THE ERROR MESSAGE

Error in sym/subs>mupadsubs (line 139)      <--- THIS FUNCTION CALLED symengine
G = mupadmex('symobj::fullsubs',F.s,X2,Y2); <--- THIS IS THE LINE THAT CALLED symengine

Error in sym/subs (line 124)                <--- THIS FUNCTION CALLED mupadsubs
G = mupadsubs(F,X,Y);                       <--- THIS IS THE LINE THAT CALLED mupadsubs

Error in nyquist2 (line 8)                  <--- THIS FUNCTION (YOURS) CALLED subs
nyq = eval(subs(G_w,w,W))                   <--- THIS IS THE LINE THAT CALLED subs


5

Suever
16 Янв 2017 в 06:22

Ответ @Suever определенно является лучшим решением в этом случае, когда необходимо вычислить большое количество значений, но другое решение, если вы хотите оценить функцию только по одному значению (0) и хотите избежать ошибки деления на ноль , ты можешь сделать это:

>> limit(G_w,w,0)
 
ans =
 
NaN

Однако это требует больших вычислительных ресурсов, поэтому его стоит использовать только тогда, когда вы ожидаете деление на ноль (например, при w = 0).


0

Bill
27 Апр 2021 в 06:24

Источник

@thilinarmtb, actually, in debug mode you cannot construct a rational that is not canonicalized. Your code:

mpq_class q = mpq_class(2, 4);
r = Rational::from_mpq(q);
std::cout << *r << std::endl;

does not construct the rational, but rather raises an exception:

test_basic: /home/ondrej/repos/csympy/src/rational.cpp:8: CSymPy::Rational::Rational(mpq_class): Assertion `is_canonical(this->i)' failed.

Abort caught. Printing stacktrace:

Traceback (most recent call last):
  File unknown, in _start()
  File "/build/buildd/eglibc-2.19/csu/libc-start.c", line 287, in __libc_start_main()
  File "/home/ondrej/repos/csympy/src/tests/basic/test_basic.cpp", line 454, in main()
    test_print_minus_one();
  File "/home/ondrej/repos/csympy/src/tests/basic/test_basic.cpp", line 429, in test_print_minus_one()
    r1 = Rational::from_mpq(q);
  File "/home/ondrej/repos/csympy/src/rational.cpp", line 29, in CSymPy::Rational::from_mpq(__gmp_expr<__mpq_struct [1], __mpq_struct [1]>)
    return rcp(new Rational(i));
  File "/home/ondrej/repos/csympy/src/rational.cpp", line 8, in CSymPy::Rational::Rational(__gmp_expr<__mpq_struct [1], __mpq_struct [1]>)
    CSYMPY_ASSERT(is_canonical(this->i))
  File "/build/buildd/eglibc-2.19/assert/assert.c", line 101, in __GI___assert_fail()
  File "/build/buildd/eglibc-2.19/assert/assert.c", line 92, in __assert_fail_base()
  File "/build/buildd/eglibc-2.19/stdlib/abort.c", line 89, in __GI_abort()
  File "/build/buildd/eglibc-2.19/signal/../nptl/sysdeps/unix/sysv/linux/raise.c", line 56, in __GI_raise()
  File unknown, in killpg()
  File "/home/ondrej/repos/csympy/src/teuchos/Teuchos_stacktrace.cpp", line 447, in loc_abort_callback_print_stack()
    Teuchos::show_stacktrace();

So internally, you don’t have to worry, because any such code would not pass tests (in debug mode) on Travis. However, you are right, that users (e.g. through Python wrappers) must not be allowed to trigger such behavior, since the production version is compiled in the Release mode.

As to @sushant-hiray’s original code, that indeed divides by zero even in debug mode, so that needs to be fixed.

Источник

Понравилась статья? Поделить с друзьями:
  • Error using save unable to write file permission denied
  • Error using plot vectors must be the same length матлаб
  • Error using plot data must be numeric datetime duration or an array convertible to double
  • Error using open line 162 matlab
  • Error using not enough input arguments matlab