Finding the roots of a high degree equation?

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)

[Torsten]

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]

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

Answer 1

[Walter Roberson]

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.

[SAM]

My problem is that I need to draw a graph where L is the x-axis so if I spacified one value for L, I won't be able to do the graph.

By smallest I mean real positive smallest value.

Anyways thank you

[Torsten]

Make a loop over all numerical L values of interest, solve your problem for each L separately and make the graph.

[Walter Roberson]

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.

[Torsten]

EI and k are given. Second last line of the OP's code.

[Walter Roberson]

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]

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]

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