# Finding the roots of a high degree equation?

1 view (last 30 days)
SAM on 20 Nov 2022
Edited: Torsten on 23 Nov 2022
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 CommentsShowHide 1 older comment
Jan on 20 Nov 2022
@SAM: This is not twitter - no # before tha tags. Thanks.

Walter Roberson on 20 Nov 2022
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.
Torsten on 23 Nov 2022
Edited: Torsten on 23 Nov 2022
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

### Categories

Find more on Function Creation in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by