I have already successfully run a code for the simulation of the deflection of beams under different loadings. I used the Matlab program and the ode45 solver for Initial value Problems and bvp4c solver for Boundary conditions. So I generated graphics such as the one below:
I want now to plot not only the deflection y(x) but also dy/dx and d^2y/dx^2 on the same graphic and I am having difficulties to actually do it. Hope someone can give me a hint. All I want is to plot the function I inserted on the MyFunctionL and MyFunctionNL and it´s derivative together with the result. One of the codes I generated is below:
function def1= def1(F,L)
%--------------------------------------------------------------
% Deflection
% F [N] und L [mm]
% EI N*mm^2
% y[mm]
%--------------------------------------------------------------
global Fg Lg EI;
Fg = F;
Lg = L;
Em=200*10^9
%This part below is not relevant for the question
for i=1:3
if i==1
b=0.055
h=0.1
Im=b*h^3/12
end
if i==2
Im=2.5175*10^(-6)
end
if i==3
re=0.065
ri=0.060
Im=pi/4*((re^4)-(ri^4))
end
%As Im is in m^4 we are converting EI to N*mm^2 by multiplying it by (1000^2).
EI=(Em*Im)*(1000^2)
$Now this part below is.
%Start point
x0 = [0 0];
%Längenintervall
xspan = [0 L];
%Call Solver -> Linear
[x y] = ode45(@MyFunctionL,xspan, x0);
%Plot result
figure
plot(x,y(:,1));
if i==1 title('y(x) in a Retangular Profile')
end
if i==2 title('y(x) in a Iprofile(IPE 100)')
end
if i==3 title('y(x) in a Round hollow section')
end
set(gca,'ydir','reverse')
xlabel('[mm]');
ylabel('[mm]');
hold on;
%Call Solver -> NonLinear
[x z] = ode45(@MyFunction,xspan, x0);
%Plot result
plot(x,z(:,1),'-r');
set(gca,'ydir','reverse')
end
return
%---------------------------------------------------------------
function dy = MyFunctionL(x,y)
global Fg Lg EI;
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = (Fg/EI)*(Lg - x);
return
function dy = MyFunction(x,y)
global Fg Lg EI;
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = (Fg/EI)*(Lg - x)*(1+y(2)^2)^1.5;
return
Best Answer
There is no need to use inaccurate finite difference approximation. You already have an exact expression for your derivatives: the functions
MyFunctionL
andMyFunction
. You simply need to evaluate them using thex
andy
values obtained from numerical integration withode45
. To do this you can either vectorize these functions or wrap them infor
loops. Additionally, you should not be using global variables. That's not what they're for and they slow down your code. You should pass additional parameter via an anonymous function. Here's how you might re-write your two integration functions to vectorize them and not use globals:You can then call the sub-functions (they could be in separate M-files too), passing in the additional parameters, and obtain the derivatives like this:
Note the transposes (
.'
) which are present becauseode45
needs a column vector but outputs each time step along the columns (for performance reasons). You could equally write separate functions that use your expressions for the derivates.