Finding the roots of a high degree equation?

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

Sign in to comment.

Accepted Answer

Walter Roberson
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.
  7 Comments
Torsten
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

Sign in to comment.

More Answers (0)

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