contour plot real numbers

I wrote a program that just takes an input of x and y values and calculates the principal stresses. It keeps giving me an error about input arguments for contour must be real. All the variables are squared so I am not sure how I would be getting this problem. Thanks.

clc

clear all

% Load

P = 1;

% Length of Beam

L = 1;

% Height of Beam from Neutral Axis

c = 0.1;

x = linspace(0,L);

y = linspace(-c,c);

[X,Y] = meshgrid(x,y);

sigmaP1 = ((3*P/(4*(c^3)))*[X*Y+[((c^2-Y^2)^2)+(X^2)*(Y^2)]^(1/2)]);

contourf(X,Y,sigmaP1,'ShowText','on')

Answer 1

[Adam Danz]

The values of 'sigmaP1 ' are imaginary. You can identify this just by looking at those values.

sigmaP1 (1:5)
ans =
  Columns 1 through 5
       1282.8 + 1.9693e-08i  1282.8 + 8.1557e-09i  1282.8 + 2.4166e-08i  1282.8 - 2.1095e-07i  1282.8 - 2.3987e-08i

Or by using isreal()

 isreal(sigmaP1)
        ans =
          logical
           0 %which means that your data contains at least 1 imaginary number

In fact, none of the elements of sigmaP1 are real numbers as confirmed by

any(imag(sigmaP1(:))==0)
ans =
  logical
   0 %no real numbers found

[Derek Reitnouer]

I did notice this when I expanded sigmaP1. I can't figure out why it is getting imaginary numbers when all the terms are squared inside the square root part of the equations.

[Adam Danz]

If your intentions are to take the square root of every element in the matrix produced by

[((c^2-Y^2)^2)+(X^2)*(Y^2)]

you must include the dot before the exponent.

sigmaP1 = ((3*P/(4*(c^3)))*[X*Y+[((c^2-Y^2)^2)+(X^2)*(Y^2)].^(1/2)]); % note the dot .^

That will produce all real numbers. However, the 'sigmalP1' matrix will contain nearly identical values which won't provide much of a contour.

Also some of your square brackets can be replaced by parentesis but I did not make that change in the line above.

[Derek Reitnouer]

Thanks for the help so far. I must not be understanding this correctly. It should produce nice countour lines I thought. I just changed the program so a set increment. I thought the sigmaP1 should give me a matrix of sigmaP1 values tied to each x and y combination. Is this correct?

clc

clear all

% Load

P = 1;

% Length of Beam

L = 1;

% Height of Beam from Neutral Axis

c = 0.1;

delta = 20;

x = [0:L/delta:L];

y = [-c:2*c/delta:c];

[X,Y] = meshgrid(x,y);

sigmaP1 = (3*P/(4*(c^3)))*(X*Y+(((c^2-Y.^2)^2)+(X.^2)*(Y.^2)).^(1/2));

contour(X,Y,sigmaP1,'ShowText','on')

[Derek Reitnouer]

I added .*

clc

clear all

% Load

P = 1;

% Length of Beam

L = 1;

% Height of Beam from Neutral Axis

c = 0.1;

delta = 20;

x = [0:L/delta:L];

y = [-c:2*c/delta:c];

[X,Y] = meshgrid(x,y);

sigmaP1 = (3*P/(4*(c^3)))*(X.*Y+(((c^2-Y.^2)^2)+(X.^2)*(Y.^2)).^(1/2));

contour(X,Y,sigmaP1,'ShowText','on')

[Adam Danz]

The code you provided above does produce contour lines (image below). Is there anything else I can help you with?

[Derek Reitnouer]

The part I'm not understanding is when I calculate the values by hand, it doesn't match what sigmaP1 is coming up with.

[Adam Danz]

Then something is wrong in your sigmaP1 equation or in the steps you're using to calculate the values by hand. I suggest running parts of the equation independently to determine if you're getting the expected values.

For example, you can 'dissect' the equation into component parts and determine if each component part makes sense. Is it the expected shape and size? Are the values as expected? Is this how you're calculating the values by hand?

% Full equation
sigmaP1 = (3*P/(4*(c^3)))*(X.*Y+(((c^2-Y.^2)^2)+(X.^2)*(Y.^2)).^(1/2));

% Run each of these sections independently to understand how the equation 
% breaks apart into components
s1 = (3*P/(4*(c^3))); 
s2 = (X.*Y+(((c^2-Y.^2)^2)+(X.^2)*(Y.^2)).^(1/2));
sigmaP1_1 = s1 * s2; 
isequal(sigmaP1, sigmaP1_1)     %check that dissection is correct (should ==1)


s1 = (3*P/(4*(c^3))); 
s2 = X.*Y; 
s3 = (((c^2-Y.^2)^2)+(X.^2)*(Y.^2)).^(1/2); 
s4 = s2 + s3; 
sigmaP1_2 = s1 * s4; 
isequal(sigmaP1, sigmaP1_2)     %check that dissection is correct


s1 = (3*P/(4*(c^3))); 
s2 = X.*Y; 
s3 = (((c^2-Y.^2)^2)+(X.^2)*(Y.^2)); 
s4 = s3.^(1/2); 
s5 = s2 + s4; 
sigmaP1_3 = s1 * s5; 
isequal(sigmaP1, sigmaP1_3)     %check that dissection is correct



s1 = (3*P/(4*(c^3))); 
s2 = X.*Y; 
s3 = ((c^2-Y.^2)^2); 
s4 = (X.^2)*(Y.^2);
s5 = s3 + s4; 
s6 = s5.^(1/2); 
s7 = s2 + s6; 
sigmaP1_4 = s1 * s7; 
isequal(sigmaP1, sigmaP1_4)     %check that dissection is correct

[Adam Danz]

Sometimes it helps to see the equation (roughly) written out.

[Derek Reitnouer]

That equation is correct how you have it written. I been tied up with work. I want to take a closer look at what you did there. Thanks for all the help so far.

[Derek Reitnouer]

When you broke it down, some of them calculate right. Others do not. s3 and s4 I checked out. I reduced the linspace down to only 5 numbers to make checking easier.

clc

clear all

% Load

P = 1;

% Length of Beam

L = 1;

% Height of Beam from Neutral Axis

c = 0.1;

delta = 4;

x = [0:L/delta:L];

y = [-c:2*c/delta:c];

[X,Y] = meshgrid(x,y);

% Full equation

sigmaP1 = (3*P/(4*(c^3)))*(X.*Y+(((c^2-Y.^2)^2)+(X.^2)*(Y.^2)).^(1/2));

% Run each of these sections independently to understand how the equation

% breaks apart into components

s1 = (3*P/(4*(c^3)));

s2 = (X.*Y+(((c^2-Y.^2)^2)+(X.^2)*(Y.^2)).^(1/2));

sigmaP1_1 = s1 * s2;

isequal(sigmaP1, sigmaP1_1) %check that dissection is correct (should ==1)

s1 = (3*P/(4*(c^3)));

s2 = X.*Y;

s3 = (((c^2-Y.^2)^2)+(X.^2)*(Y.^2)).^(1/2);

s4 = s2 + s3;

sigmaP1_2 = s1 * s4;

isequal(sigmaP1, sigmaP1_2) %check that dissection is correct

s1 = (3*P/(4*(c^3)));

s2 = X.*Y;

s3 = (((c^2-Y.^2)^2)+(X.^2)*(Y.^2));

s4 = s3.^(1/2);

s5 = s2 + s4;

sigmaP1_3 = s1 * s5;

isequal(sigmaP1, sigmaP1_3) %check that dissection is correct

s1 = (3*P/(4*(c^3)));

s2 = X.*Y;

s3 = ((c^2-Y.^2)^2);

s4 = (X.^2)*(Y.^2);

s5 = s3 + s4;

s6 = s5.^(1/2);

s7 = s2 + s6;

sigmaP1_4 = s1 * s7;

isequal(sigmaP1, sigmaP1_4) %check that dissection is correct

[Adam Danz]

This is telling. Since my break-down produces the same exact results as the full equation, that confirms that the break-down is correct. But since the break-down results to not match your expected results, that tells us that the original equation is wrong (or your expected results are wrong).

I think it's a great idea to reduce the number of inputs to ~5 so you can manage the comparison. I haven't read through your code above and I didn't see a specific follow-up question but if you get stuck in this trouble-shooting process, let me know.

[Derek Reitnouer]

Thats correct, the breakdown does match the full equation. Maybe I'm not understanding the way meshgrid works. If I go down to 3 variables, the big X equals a 3x3 matrix consisting of 0, 0.5, 1. The big Y is the same with variables -0.1,00.1. The output for sigmP1(1,1) should be calculated using the values X=0 and Y = -0.1. The equation you wrote down is the correct equations and it should return 0, not 75.