MATLAB: Changing from ‘bar’ command to ‘plot’ command.

graphplot

Double expansion chamber program:

l1=0.05;        
l3=0.05;
l5=0.05;
l2=0.200;        
l4=0.200;        
d1=0.05;
d3=d1;
d5=d1;
d2=0.25;
d4=d2;
M=0;
S1=pi/4*d1*d1;
S3=S1;
S5=S1;
S2=pi/4*d2*d2;
S4=S2;
for z=1:4500;
k=2*pi*z/343;
SP1=[cos(k*l1) 1i*sin(k*l1); 1i*sin(k*l1) cos(k*l1)]; % 2X2 matrix





TP1=[1 0; 0 S2/S1];
SP2=[cos(k*l2) 1i*sin(k*l2); 1i*sin(k*l2) cos(k*l2)];% 2X2 matrix
TP2=[1 0; 0 S3/S2];
SP3=[cos(k*l3) 1i*sin(k*l3); 1i*sin(k*l3) cos(k*l3)];% 2X2 matrix
TP3=[1 0; 0 S4/S3];
SP4=[cos(k*l4) 1i*sin(k*l4); 1i*sin(k*l4) cos(k*l4)];% 2X2 matrix
TP4=[1 0; 0 S5/S4];
SP5=[cos(k*l5) 1i*sin(k*l5); 1i*sin(k*l5) cos(k*l5)];% 2X2 matrix
TLM=SP1*TP1*SP2*TP2*SP3*TP3*SP4*TP4*SP5; % 2X2 matrix
Sum=TLM(1,1)+TLM(2,2)+TLM(2,1)+TLM(1,2); 
Mag=abs(Sum);
TL(z)=20*log10(Mag/2); 
f(z)=z;
end;
   figure(1)
   bar(f,TL,'y');
   xlabel('FREQUENCY in Hertz');
   ylabel('TRANSMISSION LOSS IN dB');
   title('Plot of TRANSMISSION LOSS Versus FREQUENCY');

I want to do this program using 'plot' command rather than 'bar' command. Its showing some errors such as >>too many input arguments.

Best Answer

You probably have a non-builtin function called 'bar.m' or 'plot.m' somewhere on your path. Execute

 >>which bar(f,TL,'y')
 >>which plot(f,TL,'y')

in the workspace of your code to find out.

Related Solutions

MATLAB: Storing output (which are matrices) of for loops

Array preallocation works for me:

l = 1;
B = [0.0250; 0.0500; 0.0750; 0.1000; 0.1250; 0.1500; 0.1750; 0.2000; 0.2250; 0.2500; 0.2750; 0.3000; 0.3250; 0.3500; 0.3750; 0.4000];
B = 0.025:0.025:0.4;
kV = nan(1,1,numel(B));
CM = nan(2,4,numel(B));
rig = 0.01;
for j = 1:numel(B)                          
    k = B(j)/(2*rig);                     
    TMP = [cos(k*l)^2,        (sin(k*l)*cos(k*l))/k, sin(k*l)*cos(k*l),    (sin(k*l)^2)/k;      
        -k*sin(k*l)*cos(k*l), cos(k*l)^2,            -k*sin(k*l)^2,        sin(k*l)*cos(k*l); 
        -sin(k*l)*cos(k*l),   -(sin(k*l)^2)/k,       cos(k*l)^2,           (sin(k*l)*cos(k*l))/k;
        k*sin(k*l)^2,         -sin(k*l)*cos(k*l),    -k*sin(k*l)*cos(k*l), cos(k*l)^2];
    TMP = TMP*TMP;
    CM(:,:,j) = [TMP(1,:);TMP(3,:)];
    kV(j) = k;
end

MATLAB: How to run the animation of double pendulum chaotic nature

Compare the results with the ODE45 integrator:

[T, Y] = ode45(F, [0, tmax], s(:, 1).');
figure;
plot(T, Y);
figure;
plot(s.')  % Your solution

The trajectories are similar, but not equal. This is a problem of the too large step width h=0.1. As soon as you reduce the step width to h=0.01 the results become nearly equal. Your Runge Kutta seems to be fine.

while t < tmax create the last time point 10.1, not 10.0. For a constant step width prefer:

t = linspace(0, 3, round(3 / 0.1) + 1);
s = zeros(4, numel(t));   % Pre-allocate!
for i = 2:numel(t)
  ...
end

Letting an array grow iteratively requires a lot of resources nd can reduce the runtime dramatically. Therefore a pre-allocation of the output is a good programming practice.

Start the loop at time point 2, because at t(1) the initial values are defined already. The wrong last time point is a small inaccuracy, but it can be avoided easily.

But you are right: I do not trust the trajectories also. If the integration is okay, the ODE must have a problem. How did you obtain it? Compare it with e.g. http://sophia.dtp.fmph.uniba.sk/~kovacik/doublePendulum.pdf. E.g. in the 4th component you have:

 ... - g*sin(s(2))*(m1+m2) / l2*(m1+m2)-m2*l2*cos(s(1)-s(2))^2

This means that only the last component of tne numerator is divided, so a parenthesis is missing:

( ... - g*sin(s(2))*(m1+m2) ) / l2*(m1+m2)-m2*l2*cos(s(1)-s(2))^2
^                           ^

For the first term of the 3rd component you have

-m2 * l1 * (s(4))^2 ...

but in the publication there is:

-m2 * l1 * (s(3))^2 ...

Writing this formula in a a large unreadable line as an anonymous function makes it hard to debug. Don't do this. Prefer a dedicated function, which allows to split the parts and check them using the debugger:

function dF = DoublePendulum(t, s, m1, m2, l1, l2, g)
s12 = s(1) - s(2);
m12 = m1 + m2;
n1 = -m2 * l1 * s(3)^2 * sin(s12) * cos(s12) + g * m2 * sin(s(2)) * cos(s12) - ...
   m2 * l2 * s(4)^2 * sin(s12) - m12 * g * sin(s(1));
d1 = l1 * m12 - m2 * l1 * cos(s12)^2;
n2 =  m2 * l2 * s(4)^2 * sin(s12) * cos(s12) + g * m12 * sin(s(1)) * cos(s12) + ...
   m12 * l1 * s(3)^2 * sin(s12) - m12 * g * sin(s(2));
d2 = l2 * m12 - m2 * l2 * cos(s12)^2;
dF = [s(3); ...
   s(4); ...
   n1 / d1, ...
   n2 / d2];
end

Now the pattern of the two equations is visible on first sight and finding typos like the s(3) or s(4) is much easier.

Then define in your code:

F = @(t,s) DoublePendulum(t, s, m1, m2, l1, l2, g);

By the way: It is a drawback to obtain the dependent X/Y coordinates befor the integeration. Integrating the formulas for the angles is much faster and more accurate. Then the results can be used easily to compute the karthesian coordinates afterwards for the visualization. For a double pendulum both will work, but if you simulate > 4 connected pendulums the computing time with dependent coordinates explodes and the accuracy degrades.

Best Answer

Related Solutions

MATLAB: Storing output (which are matrices) of for loops

MATLAB: How to run the animation of double pendulum chaotic nature

Related Question