How to plot not only the result but also the derivatives using the ode45 function?
5 views (last 30 days)
Show older comments
Daniel Clemente
on 24 Apr 2014
Commented: Star Strider
on 24 Apr 2014
Hello there, I have already successfully run a code for the simulation of the deflection of beams under different loadings. I used the ode45 for Initial value Problems and bvp4c 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
0 Comments
Accepted Answer
Star Strider
on 24 Apr 2014
Edited: Star Strider
on 24 Apr 2014
The ODE solvers return (in your situation) [y(x) y’(x)] and [z(x) z’(x)], so the first derivative is in y(:,2) and z(:,2). I didn’t run your code, but I ran a similar version of these functions with an ODE I have been working with, and since this strategy worked with it, I adapted yours to:
MyDrvsL = @(x,y) [y(:,2)'; (Fg/EI).*(Lg - x')];
MyDrvs = @(x,y) [y(:,2)'; (Fg/EI).*(Lg - x').*(1+y(:,2)'.^2).^1.5];
Given (Nx2) vectors for y and z, these should return [y’ y”]’ as a (2xN) matrix. (Note the transposition.)
2 Comments
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!