[Math] How to plot not only the result but also the derivatives of an ode using the ode45 function in Matlab

graphing-functionsMATLABordinary differential equations

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:

enter image description here

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 and MyFunction. You simply need to evaluate them using the x and y values obtained from numerical integration with ode45. To do this you can either vectorize these functions or wrap them in for 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:

function dy = MyFunctionL(x,y,Fg,Lg,EI)
dy(1,:) = y(2,:);
dy(2,:) = (Fg/EI)*(Lg-x);

function dy = MyFunction(x,y,Fg,Lg,EI)
dy(1,:) = y(2,:);
dy(2,:) = (Fg/EI)*(Lg-x).*(1+y(2,:).^2).^1.5;

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:

[x,y] = ode45(@(x,y)MyFunctionL(x,y,Fg,Lg,EI),xspan,x0);
xdot = MyFunctionL(x.',y.',Fg,Lg,EI).';
...
[x,z] = ode45(@(x,y)MyFunction(x,y,Fg,Lg,EI),xspan,x0);
zdot = MyFunctionL(x.',y.',Fg,Lg,EI).';

Note the transposes (.') which are present because ode45 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.

Related Question