Finding the roots of a high degree equation?

Question

SAM 2022-11-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1857473-finding-the-roots-of-a-high-degree-equation

编辑： Torsten 2022-11-23

Hello,

I have the following 20*20 matrix and the determinant of the matrix equal zero. The determinant of the matrix of n degree in terms of P1. Now I know that this equation will give 20 values of p1 and I want to find the smallest value p1. Could anyone please help me with this.

syms z L n k EI p c p1 
n =20;
sz_i = n;
sz_j = n;
A = zeros(sz_i, sz_j); % Preallocate your memory
for ni=1:sz_i
    for nj=1:sz_j
        yi=(1-cos((2*ni-1)*pi/2/L*z));
        yj=(1-cos((2*nj-1)*pi/2/L*z));
        yi1=diff(yi,z);
        yi2=diff(yi1,z);
        yj1=diff(yj,z);
        yj2=diff(yj1,z);
        A1(ni, nj) =int(yi2*yj2,z,0,L)*L^2/pi^2;
        A2(ni, nj) =int(k/EI*yi*yj,z,0,L)*L^2/pi^2;
        A3(ni, nj) =int(yi1*yj1,z,0,L)*p1;
        A=A1+A2+A3;
        end
end
det=simplify(det(A)==0)
det=subs(det,[EI k],[1.4*10^6 20000])
p1=solve(det,p1)

2 个评论
显示无隐藏无

Torsten 2022-11-20

编辑：Torsten 2022-11-20

And you are sure that it is possible to decide which p1 is smallest without specifying L as a numerical value ? To be honest: I doubt it.

And if you specify L as a numerical value: Use "roots" to determine p1, not "solve".

Jan 2022-11-20

@SAM: This is not twitter - no # before tha tags. Thanks.

请先登录，再进行评论。

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

Answer 1

Walter Roberson 2022-11-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1857473-finding-the-roots-of-a-high-degree-equation#answer_1106303

在 MATLAB Online 中打开

You can increase the efficiency.

syms z L n k EI p c p1 N
Pi = sym(pi);
n =20;
sz_i = n;
sz_j = n;
Y = (1-cos((2*N-1)*Pi/2/L*z));
Y1 = diff(Y, z);
Y2 = diff(Y1, z);
a12factor = L^2/Pi^2;
a3factor = p1;
yn = subs(Y, N, 1:n);
y1n = subs(Y1, N, 1:n);
y2n = subs(Y2, N, 1:n);
a1 = y2n .* y2n.' * a12factor;
a2 = yn .* yn.' * k/EI * a12factor;
a3 = y1n .* y1n.' * a3factor;
a = a1 + a2 + a3;
A = int(a, z,0,L);
D = det(A);

This is a symbolic expression in EI, k, L, p1 . You can ask MATLAB to solve() it, which will give you a root() expression with roots 1 to 20 that starts with

291703339392448328621651929229077582620568440174548205553784771502333756969121739673547250585656738532071173191070556640625*EI^20*pi^82 - 12848203493112198395889046043812884859217084764445999104*L^80*k^20 - 208389359856899544203781763568613158007234205266739200*L^80*k^20*pi + 1424916150093398371686638159155142565090345214615234107400096054751498986800574941200459804113425927420089552307128906250000*EI^20*z*pi^82 + 1381792269762211294561843388310775683848723402588160000*L^80*k^20*pi^2 + 1112350495497485111654198329296281083177130185179996821848818709250431543376954389692834277455824804598413814331054687500000*EI^20*z^2*pi^82 + 332138129688311341361837332491557801987003745461437160020708407556544889529014042088568410231691049804871850039062500000000*EI^20*z^3*pi^82 + 50648283370987298871613343472725217281132008252741289242977036400590501580538648738773426097648731595598064907812500000000*EI^20*z^4*pi^82

The coefficients get as high as 10^163, and that tells you that the solutions are going to be rather sensitive to floating point round-off.

If you had exact values for EI, k, and L, then you could subs() into the det and solve() that, getting out a series of root() with purely numeric coefficients. You could then digits(200) and vpa() . However, it is highly likely that many of the results will be complex-valued with positive and negative real components and imaginary components, so you need to define what "smallest" means over a complex set.

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

Walter Roberson 2022-11-22

在 MATLAB Online 中打开

Note that having a particular numeric value of L to concentrate on will not solve the problem that EI and k are symbolic in the expressions.

If you have numeric EI, k, and L, then you can allRoots = vpasolve(D) to get a set of 20 roots, and then you can

minRoot = min(allRoots(imag(AllRoots) == 0 & real(AllRoots) > 0));

However note that if you have numeric EI and k, and you are hoping to have symbolic L and hoping to be able to create minRoot as a formula in L, so that you could then substitute in particular numeric L and get an answer out directly... then unfortunately that is not going to work in most cases. In some cases it is possible to show that all of the 20 roots are complex-valued except for 2, and then to further prove that one particular root will always be the minimum root (or to prove that there is a particular L value that is the breakpoint and that for lower L values it is one particular root and for higher L values it is the other root). Sometimes you can do that using theory. But MATLAB is rarely going to be able to prove such a thing itself given an arbitrary degree 20 polynomial.

Walter Roberson 2022-11-22

在 MATLAB Online 中打开

You are right, @Torsten I overlooked that.

When those particular values are substituted in before-hand,

syms z L n k EI p c p1 N
Pi = sym(pi);
EI = sym(14)*10^5;
k = sym(20000);
n = 20;
sz_i = n;
sz_j = n;
Y = (1-cos((2*N-1)*Pi/2/L*z));
Y1 = diff(Y, z);
Y2 = diff(Y1, z);
a12factor = L^2/Pi^2;
a3factor = p1;
yn = subs(Y, N, 1:n);
y1n = subs(Y1, N, 1:n);
y2n = subs(Y2, N, 1:n);
a1 = y2n .* y2n.' * a12factor;
a2 = yn .* yn.' * k/EI * a12factor;
a3 = y1n .* y1n.' * a3factor;
a = a1 + a2 + a3;
tic
A = int(a, z,0,L);
toc    %takes over 6 minutes on my machine
tic
D = det(A);
toc   %takes about 45 seconds on my machine
sol = solve(D, p1);
CH = children(sol(1),1);

then CH will be a polynomial of degree 20 in z, whos 20 roots are the potential solutions.

If you

[P,Q] = coeffs(CH,L)

you will see that all of the powers of L are 4th powers, such as L^20 L^16 . That tells you that there is no need to investigate negative L: if there are no solutions in positive L then there are no solutions in negative L either.

If you continue with

[P, Q] = coeffs(CH,z);
vpa(P)

you will see that P consists entirely of positive coefficients together with positive coefficients times L to even powers. None of the terms can potentially be negative. And that means that none of the roots can possibly be positive.

Which is to say, at least for those values of EI and k, given real-valued L, there are no possible positive roots p1 for the determinant. All of the roots turn out to be real-valued and negative. For small enough absolute value of L, the roots are

-0.25

-2.25

-6.25

-12.25

-20.25

-30.25

-42.25

-56.25

-72.25

-90.25

-110.25

-132.25

-156.25

-182.25

-210.25

-240.25

-272.25

-306.25

-342.25

-380.25

to within round-off error. Which are differences of 38, 34, 32, ... 2. As L increases, the roots get further from that pattern.

SAM 2022-11-23

在 MATLAB Online 中打开

yes, I realized that A3 should be negative. I also tried to normalize the problem which kept the determinant at the end as a factor of p1 and a. It all worked fine, but when I tried to use another for lope between a and p1 to get the different values of p1, I faced the following problem. I used the subs function to substitute r instead of a but I realized that Matlab is actually substituting the same value of r every time. for example when I let r=0:4 thats what I got. so is there another way I can do this loop?

P1 =

[0.2500, 0.2500, 0.2500, 0.2500, 0.2500]

clear;clc
syms z L n k EI p c p1 a 
n =20;
sz_i = n;
sz_j = n;
A = zeros(sz_i, sz_j); % Preallocate your memory
for ni=1:sz_i
    for nj=1:sz_j
        yi=(1-cos((2*ni-1)*pi/2/L*z));
        yj=(1-cos((2*nj-1)*pi/2/L*z));
        yi1=diff(yi,z);
        yi2=diff(yi1,z);
        yj1=diff(yj,z);
        yj2=diff(yj1,z);
        A1(ni, nj) =int(yi2*yj2,z,0,L)*L^3/pi^4;
        A2(ni, nj) =int(yi*yj,z,0,L)/L*a;
        %a=kL^4/pi^4/EI
        A3(ni, nj) =-int(yi1*yj1,z,0,L)*p1*L/pi^2;
        %p1=pL^2/pi^2/EI
        A=A1+A2+A3;
    end
end
det=det(A)==0;
ca=1;
    for r=0:.001:20
det=subs(det,a,r);
 allRoots =vpasolve(det,p1);
P1(ca)=min(allRoots);
ca=ca+1;
    end

Torsten 2022-11-23

编辑：Torsten 2022-11-23

在 MATLAB Online 中打开

clear;clc
syms z L n k EI p c p1 a 
n = 20;
sz_i = n;
sz_j = n;
A = zeros(sz_i, sz_j); % Preallocate your memory
for ni=1:sz_i
    for nj=1:sz_j
        yi=(1-cos((2*ni-1)*pi/2/L*z));
        yj=(1-cos((2*nj-1)*pi/2/L*z));
        yi1=diff(yi,z);
        yi2=diff(yi1,z);
        yj1=diff(yj,z);
        yj2=diff(yj1,z);
        A1(ni, nj) =int(yi2*yj2,z,0,L)*L^3/pi^4;
        A2(ni, nj) =int(yi*yj,z,0,L)/L*a;
        %a=kL^4/pi^4/EI
        A3(ni, nj) =-int(yi1*yj1,z,0,L)*p1*L/pi^2;
        %p1=pL^2/pi^2/EI
        A=A1+A2+A3;
    end
end
det=det(A)==0;
ca=1;
for r=0:0.001:20
    detnum=subs(det,a,r);
    allRoots =vpa(solve(detnum));
    P1(ca) = min(allRoots(abs(imag(allRoots))<1e-4));
    ca=ca+1;
end

请先登录，再进行评论。

Finding the roots of a high degree equation?

2 个评论
显示无隐藏无

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

Finding the roots of a high degree equation?

2 个评论 显示 无隐藏 无

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

2 个评论
显示无隐藏无

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