Help with a tiny project ( building a graph)

I will try again, if anyone can help I would be very grateful...

I was wondering if someone could help me with a problem I have been given. It is about deflection of beams in a course in college called numerical analysis.

Given: 𝑞0 = 1𝑁𝑚𝑚, E = 140MPa, L = 1000mm, v = -0.7mm.

Draw a graph of x in front of I for 100∙ 10^6 mm^4 ≤ I ≤ 450∙10^6 mm^4. What conclusion can .be physically deduced from the graph?

Now, here is the code for x, I found it numerically (function beam2).I used the newton raphson method.

function X_coordinate = beam2(x0,q0,v,L,I,E)

epsilon = 1e-4;
flag = 0;
counter = 0;
% v=-0.5
% x0=500;
% q0=0.5
% L=1000
% I=351*10^6
% E= 100
x = x0;
while flag == 0
    %Check for infinite loop
    counter = counter + 1
    if counter > 10000;
        error('Too many iterations');
    end
    F = -(q0*L)/(3*pi^4*E*I)*(48*L^3*cos((pi*x)/(2*L))-48*L^3+3*pi^3*L*x^2-pi^3*x^3)-v
    F_d = -(q0*L)/(3*pi^4*E*I)*(-48*L^3*(pi/(2*L))*sin((pi*x)/(2*L))+6*pi^3*L*x-3*pi^3*x^2)
     x = x-F/F_d %solution 
    
    if abs(F/F_d) < epsilon || abs(F) < epsilon 
        flag = 1;
    end
end
X_coordinate = x;
if X_coordinate <=0 || X_coordinate > 1000
    error('No solution found,Please enter valid value')
else
    X_coordinate = x
end

The code down is the code for buliding the graph. I need to use the solution of the beam2 function for that, numerically of course, and I have no idea how to proceed from here.

%function x_I = MomentI(x0,q0,L,v,E)

epsilon = 1e-4;
flag = 0;
counter = 0;
 v=-0.7
 x0=500;
 q0=1
 L=1000
% I=351*10^6
 E= 140
 x = x0
% I_x = 1:10
% x_x = 1:10

% while flag == 0
    %Check for infinite loop
%     counter = counter + 1
%     if counter > 10;
%         error('Too many iterations');
%     end
    if I_x >= 100*10^6 & I_x<=450*10^6
         for x=1:1000
            I_x(x) = -(q0*L)/(3*pi.^4*E*v)*(48*L^3*cos((pi*x)/(2*L))-48*L.^3+3*pi.^3*L*x.^2-pi^3*x.^3)
            I_d(x) = -(q0*L)/(3*pi.^4*E*v)*(-48*L.^3*(pi/(2*L))*sin((pi*x)./(2*L))+6*pi^3*L*x-3*pi.^3*x.^2)
            x_x(x) = x- I_x(x)./I_d(x) %solution 
         end
    end
%     if abs(I_x(x)./I_d(x)) < epsilon | abs(I_x(x)) < epsilon 
%         flag = 1;
%     end
    %end
% end
x_I = x;
if x_I <=0
    error('No solution found,Please enter valid value')
else
    x_I = x
end
x = 1:1000;
hold on
title ('I as a function of x')
xlabel('x[mm]')
ylabel('I[mm]')
plot(x_x,I_x)
hold off

Any suggestions to improve the current code would be great . thanks ...

Answer 1

[Chris]

If I understand correctly, you only wish to plot I as a function of x.

I have modified your function slightly (at the end of the code), including some semicolons to make it less talkative.

% Make a vector of I values
I_x = linspace(100e6,450e6,1000);

% Set constants
v=-0.7;
x0=500;
q0=1;
L=1000;
E=140;

% Preallocate x vector for speed
x_x = zeros(size(I_x));

% Repeatedly call the function, assigning the result to an element of the
% x_x vector
for idx = 1:numel(I_x)
    x_x(idx) = beam2(x0,q0,v,L, I_x(idx) ,E);
end

% Plot
hold on
title ('I as a function of x')
xlabel('x[mm]')
ylabel('I[mm]')
plot(x_x,I_x)
hold off

% Your function
function X_coordinate = beam2(x0,q0,v,L,I,E)

epsilon = 1e-4;
flag = 0;
counter = 0;

x=x0;
while flag == 0
    %Check for infinite loop
    counter = counter + 1;
    if counter > 10000
        error('Too many iterations');
    end
    F = -(q0*L)/(3*pi^4*E*I)*(48*L^3*cos((pi*x)/(2*L))-48*L^3+3*pi^3*L*x^2-pi^3*x^3)-v;
    F_d = -(q0*L)/(3*pi^4*E*I)*(-48*L^3*(pi/(2*L))*sin((pi*x)/(2*L))+6*pi^3*L*x-3*pi^3*x^2);
     x = x-F/F_d; %solution 
    
    if abs(F/F_d) < epsilon || abs(F) < epsilon 
        flag = 1;
    end
end
X_coordinate = x;
if X_coordinate <=0 || X_coordinate > 1000
    error('No solution found,Please enter valid value')
else
    X_coordinate = x;
end

end