Root finding and plotting
显示 更早的评论
How do I find the three roots of " f "
I got the plot of " f." And because the " f " and X asis has three intersections, it would have three roots.
I would like to know how to find the three roots at the same time and plot the three roots on my original plot.
Below is my coding.
Thank you very much!!
clc
clear
format long
v=0.3;
E=209e+3;
G=E/(2*(1+v));
q=-1;
h=15;
D=(E*h^3)/(12*(1-v^2));
I=(h^3)/12;
a=600;b=2400;
n =3;
[T1, T2] = meshgrid(1:2:n);
mn = [T1(:), T2(:)]
syms x y
x_value=301;
y_value=0:1:2400;
len=length(mn);
amn=zeros(1,len);
for i=1:len
m=mn(i,1);
n=mn(i,2);
amn(i)=(16*q/(m*n*D*pi^6))*(1/((m/a)^2+(n/b)^2)^2);
end
wmn=sym(zeros(1,len));
wxx=sym(zeros(1,len));
wyy=sym(zeros(1,len));
for i=1:len
m=mn(i,1);
n=mn(i,2);
wmn(i)=amn(i).*((sin(m.*pi.*x./a)).*(sin(n.*pi.*y./b)));
wxx(i)=diff(wmn(i),x,2);
wyy(i)=diff(wmn(i),y,2);
end
combine_wxx=sum(wxx);
combine_wyy=sum(wyy);
My=-D*(combine_wyy+v*combine_wxx);
Differentiation_My=diff(My,y,1)
f=subs(Differentiation_My,x,301)
%fsolve(f,500)
Differentiation_My_value=double(subs(Differentiation_My,{x,y},{x_value,y_value}));
plot(Differentiation_My_value)
grid on
采纳的回答
One way is to search for the indices nearest the zero-crossings, then interpolate (if necessary, to get more exact values) —
format long g
v=0.3;
E=209e+3;
G=E/(2*(1+v));
q=-1;
h=15;
D=(E*h^3)/(12*(1-v^2));
I=(h^3)/12;
a=600;b=2400;
n =3;
[T1, T2] = meshgrid(1:2:n);
mn = [T1(:), T2(:)]
syms x y
x_value=301;
y_value=0:1:2400;
len=length(mn);
amn=zeros(1,len);
for i=1:len
m=mn(i,1);
n=mn(i,2);
amn(i)=(16*q/(m*n*D*pi^6))*(1/((m/a)^2+(n/b)^2)^2);
end
wmn=sym(zeros(1,len));
wxx=sym(zeros(1,len));
wyy=sym(zeros(1,len));
for i=1:len
m=mn(i,1);
n=mn(i,2);
wmn(i)=amn(i).*((sin(m.*pi.*x./a)).*(sin(n.*pi.*y./b)));
wxx(i)=diff(wmn(i),x,2);
wyy(i)=diff(wmn(i),y,2);
end
combine_wxx=sum(wxx);
combine_wyy=sum(wyy);
My=-D*(combine_wyy+v*combine_wxx);
Differentiation_My=diff(My,y,1)
f=subs(Differentiation_My,x,301)
%fsolve(f,500)
Differentiation_My_value=double(subs(Differentiation_My,{x,y},{x_value,y_value}));
zxi = find(diff(sign(Differentiation_My_value))) % Zero-Crossing Indices
for k = 1:numel(zxi)
idxrng = max(zxi(k)-2,1):min(zxi(k)+2,numel(Differentiation_My_value)); % Index Range For Interpolation
zc(k) = interp1(Differentiation_My_value(idxrng),idxrng,0); % Interpolate
end
Zero_Crossings = zc
plot(Differentiation_My_value)
hold on
plot(zc, zeros(size(zc)), 'rs')
hold off
grid on
To plot it against an ‘x’ vector instead of the indices, the ‘zc’ calculation becomes:
% zc(k) = interp1(Differentiation_My_value(idxrng), x(idxrng), 0);
I commented it here because it would otherwise execute in the online Run feature, and throw an error.
.
8 个评论
Mark
2021-6-17
Star Strider
2021-6-17
As always, my pleasure!
One of them is a near-duplicate (1200, 1201). I have no idea how that arises unless there is ‘noise’ or a duplicated value in that region that causes a second sign reversal. (I did not look closely at the data.) Using the unique (or uniquetol) function would remove it.
Mark
2021-6-17
Star Strider
2021-6-17
Sure!
This assignment:
idxrng = max(zxi(k)-2,1):min(zxi(k)+2,numel(Differentiation_My_value));
sets the beginning of the index range to the greater of 1 or the index value - 2, so that it does not address indices less than 1. The second part does the same at the opposite end, taking the lesser of the index value + 2 or the number of elements in the vector, so that it does not address indices greater than the number of elements in the vector. It is simply ‘crash-proofing’ the indexing range. (It is not strictly necessasry in this instance, however experience has taught me to always include it.)
‘How should I fix it if I Want to plot it against an ‘x’ vector instead of the indices.’
Since ‘x’ would have to be the same length as ‘Differentiation_My_value’, one way would be:
xstart = 0;
xend = 1;
x = linspace(xstart, xend, numel(Differentiation_My_value));
with ‘xstart’ and ‘xend’ being any values you want (within the ability of MATLAB to represent them).
As always, my pleasure!
.
Mark
2021-6-18
Star Strider
2021-6-18
As always, my pleasure!
If you have any further questions or concerns, please post them. I’ll do my best ro respond.
Mark
2021-6-18
Star Strider
2021-6-18
I’ve definitely been there (although by now a few decades removed).
As always, my pleasure!
更多回答（1 个）
Sulaymon Eshkabilov
2021-6-17
There are a couple of ways to find roots of this eqn:
(1) ginput() --> graphical method, e.g.:
[Roots, y] = ginput(3) % click on three crossing points
(2) fsolve()
类别
在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
