**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Coding an Equation to be solved

2 views (last 30 days)

Show older comments

### Accepted Answer

Alan Stevens
on 7 Apr 2021

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 Comments

Arjun Siddharth
on 7 Apr 2021

Thank you so much! This was very helpful!

Arjun Siddharth
on 7 Apr 2021

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
on 7 Apr 2021

Just dphiadx(L)

Alan Stevens
on 7 Apr 2021

It should have worked! What was the error message?

Alan Stevens
on 7 Apr 2021

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
on 7 Apr 2021

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
on 7 Apr 2021

Arjun Siddharth
on 7 Apr 2021

Alan Stevens
on 7 Apr 2021

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
on 7 Apr 2021

thanks so much!

Alan Stevens
on 8 Apr 2021

No, sorry! I must have been half-asleep! Use

phiasq = @(x) phia(x).^2;

Arjun Siddharth
on 8 Apr 2021

that's alright! thanks a lot !

Arjun Siddharth
on 8 Apr 2021

Alan Stevens
on 8 Apr 2021

Alan Stevens
on 8 Apr 2021

Edited: Alan Stevens
on 8 Apr 2021

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
on 8 Apr 2021

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
on 9 Apr 2021

Arjun Siddharth
on 9 Apr 2021

oh right. What might be causing this variation? Some other constant?

Alan Stevens
on 9 Apr 2021

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.

Arjun Siddharth
on 9 Apr 2021

alright, thanks a bunch!

Alan Stevens
on 14 Apr 2021

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
on 20 Apr 2021

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
on 20 Apr 2021

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
on 20 Apr 2021

Edited: Alan Stevens
on 20 Apr 2021

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
on 20 Apr 2021

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
on 21 Apr 2021

Alan Stevens
on 21 Apr 2021

Arjun Siddharth
on 21 Apr 2021

oh okay. thanks!

### More Answers (0)

### See Also

### Categories

### Tags

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- 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)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)