現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
Coding an Equation to be solved
1 回表示 (過去 30 日間)
古いコメントを表示
I would like some help in coding equation 14, in order to solve equation 15 to obtain the value of C. Is this possible to do? I am not able to do it.
採用された回答
Alan Stevens
2021 年 4 月 7 日
This should do it, assuming you know the values of the other constants:
% Arbitrary data
L = 1;
rho = 1;
A = 1;
sigma = 1;
Mt = 1;
It = 1;
lambda = 1.162;
k = lambda/L;
% define functions without C
phia = @(x) (cos(k*x)-cosh(k*x)+sigma*(sin(k*x)-sinh(k*x)));
phiasq = @(x) phia(x);
dphidxa = @(x) -k*(sin(k*x)+sinh(k*x)-sigma*(cos(k*x)-cosh(k*x)));
S = rho*A*integral(phiasq,0,L) + Mt*phiasq(L) + It*dphidxa(L)^2;
C = sqrt(1/S);
disp(C)
29 件のコメント
Arjun Siddharth
2021 年 4 月 7 日
Can you also tell me how to evaluate a function at a point?
For example, how do i substitute x=L in dphidxa equation?
Alan Stevens
2021 年 4 月 7 日
You just need to define beta_alpha1 as
beta_alpha1 = dphidxa(L);
It is just a number - don't make it a function!
Arjun Siddharth
2021 年 4 月 7 日
thanks! it worked!
in the intial code, could you clarify why you used 'C=sqrt(1/S)'? Is that from the mathematical solving?
Alan Stevens
2021 年 4 月 7 日
Yes. each of the terms in eqn (15) is implicitly multiplied by C^2. Hence it could be factored out. You must therefore remember to multiply the "bare" functions, phia and dphidxa, by C when you use them.
Arjun Siddharth
2021 年 4 月 7 日
how do we multiply a value like 'C' to a function? because it says operators like '*' are incompatible with function handles.
Alan Stevens
2021 年 4 月 7 日
Example:
a = 1:5;
b = C*phia(a); % not b = C*phia
You multiply C by evaluated values of the function, not the function definition itself.
Arjun Siddharth
2021 年 4 月 8 日
Hey Alan, would you have time to take a look at my code? I'm trying to recreate the graph from a research paper. The code is running, but the graph does not match. Would appreciate if you could have a look. Do let me know , so that I may send you the paper as well as the code.
Alan Stevens
2021 年 4 月 8 日
Upload your code here. Make use of the > button in the CODE section in the menu header rather than posting a picture.
Alan Stevens
2021 年 4 月 8 日
編集済み: Alan Stevens
2021 年 4 月 8 日
Your code has frequency in Hz, whereas the equations want omega in radians/sec. Your Voltage Response section is better expressed as
%% Voltage Response
freq = linspace(20,140,1200)*2*pi;
v = zeros(numel(freq),1);
for count = 1:numel(freq)
w = freq(count);
t1=-1j*w*alpha*RL*M;
t2=M*(1j*2*Deq*wn*w+wn^2-w^2);
t3=1j*w*Cp*RL+1;
t4=1j*w*RL*alpha^2;
v1=t1/(t2*t3+t4);
v(count)=abs(v1);
end
figure(1)
plot(freq/(2*pi),v,'b','Linewidth',2);
xlabel('Frequency (Hz)');
ylabel('Voltage');
The magnitude of the voltage response is still different. This might be due to differences in data - I haven't checked everything!
Your value for 'a' is 0.8558 which doesn't match the values used in figure 6.
Alan Stevens
2021 年 4 月 8 日
Your code needed some reordering and addition of some arguments in function calls. However, end result seems to be the same!
%Code to validate Improved Lumped Parameter Model
%% Parameters
rp=7800; %Density of Piezo Plate
rs=9000; %Density of Substrate Plate
Ep=66*10^9;
Es=105*10^9;
d31=-190*10^-12;
e31=-11.5;
e0=8.854*10^-12;
e33=1500*e0;
L=50.8*10^-3;
b=31.8*10^-3;
hp=0.26*10^-3;
hs=0.14*10^-3;
DrB=0.002;
Deq=0.027;
Mt=0.012;
RL=1000;
g = 9.81; % You hadn't defined g (though result seems independent of g).
%% Computed Parameters
I1=2*b*(((hp+(hs/2))^3)-((hs/2)^3));
I2=I1/3;
I3=(b*(hs^3))/12;
I=I2+I3;
E1=Es*((b*(hs^3))/12);
E2=Ep*2*b*(((hp+(hs/2))^3)-((hs/2)^3));
E3=E2/3;
E=(E1+E3)/I;
Cp=(b*L*e33)/(2*hp);
rho=((2*rp*hp)+(rs*hs))/((2*hp)+hs);
A=b*(hs+(2*hp));
It=Mt* (L^2);
%% Bending Mode Shape
lambda = 1.162;
sigma=sigma_calc(lambda,Mt);
k = lambda/L;
%Define functions without 'C'
phia = @(x) ((cos(k*x)-cosh(k*x))+(sigma*(sin(k*x)-sinh(k*x))));
phiasq = @(x) phia(x).^2;
dphidxa = @(x) -k*( sin(k*x)+sinh(k*x) - sigma*(cos(k*x)-cosh(k*x)) );
S = (rho*A*(integral(phiasq,0,L))) + (Mt*phia(L).^2) + (It*((dphidxa(L)).^2));
C = sqrt((1/S));
beta_alpha=(-L)*C*dphidxa(L);
disp(beta_alpha);
disp(C);
alpha=(beta_alpha*b*(hs+hp)*e31)/(2*L);
%% Correction Factors Computation
mew_x=@(x) (Mt*g*x.^2.*(3*L-x))./(6*E*I*C*phia(x)); % need to call phia with x as argument
mew_barx=@(x) ((-2.7*10^-5)/(99*L))*(x-(0.01*L))-(1.207*10^-5);
func= @(x)mew_barx(x).*phia(x)*C; % need to call mew_bar and phia with x as argument
func_sq=@(x) func(x).^2;
BM1=integral(func_sq,0,L);
BM2 = func_sq(L)*L;
Beta_M=BM1/BM2;
BK1=L^3*Mt*g;
BK2=E*I*mew_x(L)*C*phia(L);
Beta_K=BK1/BK2;
M_beam=(Beta_M*rho*A*L);
M=(Beta_M*rho*A*L)+Mt;
K=(Beta_K*E*I)/(L^3);
wn=sqrt((K/M));
a=Mt/M;
%% Voltage Response
freq = linspace(20,140,1200)*2*pi;
v = zeros(numel(freq),1);
for count = 1:numel(freq)
w = freq(count);
t1=-1j*w*alpha*RL*M;
t2=M*(1j*2*Deq*wn*w+wn^2-w^2);
t3=1j*w*Cp*RL+1;
t4=1j*w*RL*alpha^2;
v1=t1/(t2*t3+t4);
v(count)=abs(v1);
end
figure(1)
plot(freq/(2*pi),v,'b','Linewidth',2);
xlabel('Frequency (Hz)');
ylabel('Voltage');
%% Function to calculate Sigma Value
function [sig]=sigma_calc(lambda,Mt)
mL=0.02;
t2=lambda*Mt;
n1=sin(lambda)-sinh(lambda);
n2=cos(lambda)-cosh(lambda);
d1=cos(lambda)+cosh(lambda);
d2=sin(lambda)-sinh(lambda);
b1=(mL*n1)+(t2*n2);
b2=(mL*d1)-(t2*d2);
sig=b1/b2;
end
Alan Stevens
2021 年 4 月 9 日
You can get a = 10 by setting Mt = 0.0165. However, it only leads to a small increase in the peak; nothing like enough to match that in the Wang paper.
Alan Stevens
2021 年 4 月 9 日
Possibly some other constant.
Possibly, the graphs in the paper were generated with an entirely different set of constants, or there is an error in one or more of their published equations or their non-published coding. Since the authors don't list their coding it's not possible to know!
Perhaps there are other papers you can check against.
Alan Stevens
2021 年 4 月 14 日
I've found one error: your definition of zbar is missing some brackets. I think it should be
z_bar=(ts^2+tm^2*n+2*ts*tm*n)/(2*(ts+(n*tm)));
As a general comment, you use too many brackets! It makes the code difficult to read. You don't need them around parameters that are multiplied together, for example.
Also, you need to include more comments, explaining what the statements are, or where they connect with the paper - especially if you want help from others! I followed it down to your computed parameters section, but lost the connection with the paper at that point.
Alan Stevens
2021 年 4 月 20 日
You can fit a straight line as follows:
L = 0.01:0.01:0.05;
mu = [-1.2121, -1.2182, -1.2235, -1.2292, -1.2338];
mc0 = [-1, -1];
mc = fminsearch(@(mc) fn(mc,L,mu), mc0);
m = mc(1); c = mc(2);
disp([m,c])
-0.5440 -1.2070
x = 0.01:0.001:0.05;
mufit = m*x + c;
plot(L,mu,'o',x,mufit),grid
xlabel('x'), ylabel('\mu')
legend('data','fit')
text(0.015,-1.227,['\mu = ' num2str(m) 'x ' num2str(c)])
function F = fn(mc, L, mu)
m = mc(1); c = mc(2);
F = norm(m*L + c - mu);
end
I'll let you manipulate it to look like the expression in the JPG if that's what you want it to look like.
Alan Stevens
2021 年 4 月 20 日
The first two lines of the listing above are simply the data taken from Table 1 of the Wang paper. The resulting straight line fit for mu is mufit = m*x + c, where m and c are given as -0.54398 and -1.207 respectively, and x represents the values of length. You could define the straight line as a function:
mufit = @(x) -0.54398*x - 1.207;
then call it with the required value of x when you need to, e.g.
y = mufit(0.025);
Alan Stevens
2021 年 4 月 20 日
編集済み: Alan Stevens
2021 年 4 月 20 日
You need to recalculate a few values of mu (using equation (16)) over the range of interest of values of L, then do a new curve-fit with those values.
Note: I left out the 10^-5 from the fitted values of mu. The final fit should be multiplied by this.
Alan Stevens
2021 年 4 月 20 日
x is a dummy variable and represents values if L.
You can see from the graph that the curve fit is pretty good. It doesn't necessarily match exactly at the tabulated points because it is a best-fit line - it's not forced to go through every point. When you say the value is different, what sort of doifference do you get?
Incidentally, I should really have used the function 'polyfit' to get the line; however, as you can see below this doesn't make a significant difference here.
L = 0.01:0.01:0.05;
mu = [-1.2121, -1.2182, -1.2235, -1.2292, -1.2338];
p = polyfit(L,mu,1);
m = p(1); c = p(2);
disp([m,c])
-0.5440 -1.2070
x = 0.01:0.001:0.05;
mufit = m*x + c;
plot(L,mu,'o',x,mufit),grid
xlabel('L'), ylabel('\mu')
legend('data','fit')
text(0.015,-1.227,['\mu = ' num2str(m) 'x ' num2str(c)])
What I've called 'mufit' is the average value function; your 'mew_bar'. (Again, I've left out the 10^-5 multiplier, which is needed when you actualy use mufit.)
Alan Stevens
2021 年 4 月 21 日
If you change EI and L you should really recalculate equation 16 for a range of values of L, covering the range of interest. Here, you might manage by simply rescaling the average value of mu function by (EIold/EInew), since mu is inversely proportional to EI.
Alan Stevens
2021 年 4 月 21 日
Well, strictly, your value of L is outside the range of the values used for the curve fit, but the line is so straight I suspect this won't introduce a significant error here.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Spectral Measurements についてさらに検索
タグ
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の MathWorks のサイトは、お客様の地域からのアクセスが最適化されていません。
南北アメリカ
- América Latina (Español)
- Canada (English)
- United States (English)
ヨーロッパ
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
アジア太平洋地域
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)