How to plot not only the result but also the derivatives using the ode45 function?

5 visualizzazioni (ultimi 30 giorni)
Daniel Clemente
Daniel Clemente il 24 Apr 2014
Commentato: Star Strider il 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

Accedi per rispondere a questa domanda.

Risposta accettata

Star Strider
Star Strider il 24 Apr 2014
Modificato: Star Strider il 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.)

Più risposte (0)

Categorie

Scopri di più su Ordinary Differential Equations in Help Center e File Exchange

Tag

Prodotti

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by