Someone that works with symbolic variables?
조회 수: 2 (최근 30일)
이전 댓글 표시
Hi guys, I have a code that should give me some results. Unfortunately when I insert k>3, the program collapses and says there isn't enought RAM to run it. I'm just texting to know if anyone has ever worked with symbolic variables before and can help me. If anyone knows how to work with this please I'll be really glad to talk to you via email or Telegram.
function [S,Gnew,Enew] = BGroebner1(k)
M = 3*k + 1;
x = sym('x');
r = sym('r');
c = sym('c', [1 M+1]);
% Define the sum of monomials with negative coefficients
CC = x;
for i = 1:M+1
CC = CC + c(1,i)*x^(-i);
end
sM = expand(CC^(3));
for i = M:1:3*M+3
sM = subs(sM, x^(-i), r);
end
sMnew = sM;
sMnew = subs(sMnew,r,0);
sMnew1 = expand(x^(M-1)*sMnew);
C =coeffs(sMnew1,x);
E =sym('E',[1,M-1]);
for i= M-1:-1:1
E(1,M-i) = C(1,i);
end
% Calculate Groebner basis,
num_syms = M+1;
syms_cell = cell(1, num_syms);
for i = M+1:-1:1
syms_cell{i} = sym(['c', num2str(M+2-i)]);
end
c = [syms_cell{:}];
G = gbasis(E,c,"MonomialOrder",'Lexicographic');
J = (M+1)*M;
sM1 = expand(CC^{M});
sMnew2 = expand(x^(J)*sM1);
E2 = coeffs(sMnew2,x);
EM = E2(J);
EM1=E2(J-1);
F= sym('F');
Enew = [E,EM,EM1+F];
% Groebner basis with (E^M)_1, (E^M)_2 + F
Gnew = gbasis(Enew,[F c],"MonomialOrder",'Lexicographic');
%% Calculate the solutions of the E_{1}, E_{2},...
eqn =sym('eqn', [1 M-1]);
for i=1:1:M-1
eqn(i)= Enew(i)==0;
end
d = sym('c', [1 M+1]);
S = vpasolve(eqn,d);
end
Looking forward your replies.
댓글 수: 0
답변 (2개)
John D'Errico
2024년 6월 18일
편집: John D'Errico
2024년 6월 18일
When k is even slightly large (apparently 4 is large enough) things go to hell. Did I not answer this question before? I know I wrote the answer, since I computed eqn when k=2. In fact, I even set my computer running on the case when k==4. I broke out of the computations before it melted into a pile of hot slag though.
eqn =
[3*c1^2 + 3*c3 == 0, 3*c4 + 6*c1*c2 == 0, c1^3 + 6*c3*c1 + 3*c2^2 + 3*c5 == 0, 3*c2*c1^2 + 6*c4*c1 + 3*c6 + 6*c2*c3 == 0, 3*c1^2*c3 + 3*c1*c2^2 + 6*c5*c1 + 6*c4*c2 + 3*c3^2 + 3*c7 == 0, 3*c4*c1^2 + 6*c3*c1*c2 + 6*c6*c1 + c2^3 + 6*c5*c2 + 3*c8 + 6*c3*c4 == 0]
That is only for k=2. When k is larger, we see a far more messy, far less managable system of nonlinear polynomial equations. Any computations with those equations for values of k as large as k=4 prbably have many thousands, if not millions of terms in them.
Now, do you have an infinitely large, infinitely fast computer? Probably not, at least, not unless you are a character on a TV show or in the movies, where computers always are immensely powerful.
It is very easy to write code for an insolvable problem, or even one that is practically insoluble today on any computer you will find in your lap or a desktop. I'm sorry, but there is no magic to be found here. (And please don't send me e-mail pleading for direct e-mail consiulting. I'm saying this because you seem to be asking for someone to perform that service for you.)
What should you do? Well, you might turn the problem into one that fsolve can handle. vpasolve still uses symbolic computations, and they wil get incredibly messy, and incredibly CPU intensive. But since fsolve will be working purely in double precision arithmetic, it will have at least a chance of finding a solution. The problem is, the system you are computing will still be immensely messy. Just generating the equations is itself a computational nightmare. So you will not want to perform all of those computations using symbolic tools. You will need to work in floating point all the way. That means you need to completely reformulate how this problem will be handled.
I'm sorry, but the easy problems were all solved many years ago.
Walter Roberson
2024년 6월 18일
When k = 4,
size(char(sM1))
ans =
1 945926142
and it takes numerous minutes just to calculate the char(sM1) .
You have an explosion in terms that makes calculations impractical.
G = gbasis(E,c,"MonomialOrder",'Lexicographic');
The results in G are not used, so there is no point in computing that.
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Eigenvalue Problems에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!