MATLAB: Ode45 and diff function problem

differential equationsode45

Hello.
I made the fallowing program for project.It's supposed to show the sigma graphic with the help of ode45 function.The problem is i have a bit of a complex system whitch i need to solve.Something like:
p2*d^2sigma/dt^2+p1*dsigma/dt +p0sigma=q2*d^2epsilon/dt^2+q1*depsilon/dt+q0*epsilon (differential reprezentation of series pairing of Maxwell and Kelvin reologic models)
I wanted to make it simple at first and took the epsilon function as a single constant variable but i have some errors in useing the diff function in the "sistem_cu_ec_de_gradul_1"function.
This is the fallowing program:
function serie_maxwell_kelvin
E1=input('E1=');
niu1=input('niu1=');
E2=input('E2=');
niu2=input('niu2=');
   p2=niu2/E1;
  p1=1+E2/E1+niu2/niu1;
  p0=E2/niu1;
  q2=niu2;
  q1=E2;
  q0=0;
   t=linspace(0,10,100);
   sigma_initial=0;
   dsigmadt_initial=0;
   [t,x]=ode45(@sistemul_cu_ec_de_gradul_1,t,[sigma_initial dsigmadt_initial]);
   plot(t,x(:,1));
   xlabel('t');
   ylabel('sigma');
   function[dxdt]=sistemul_cu_ec_de_gradul_1(t,x)
     dxdt_2=x(1);
     dxdt_1=-(p1/p2)*x(1)-(p0/p2)*x(2)+(q2/p2)*diff(epsilon,t,2)+(q1/p2)*diff(epsilon,t)+(q0/p2)*epsilon(t);
     %dxdt_1=-(p1/p2)*x(1)-(p0/p2)*x(2)+(q2/p2)*0+(q1/p2)*0+(q0/p2)*epsilon(t);
     dxdt=[dxdt_2; dxdt_1];
   end
   function [z] = epsilon(t)
   z=6;
   end
   end
As start i took the E1=10,niu1=5,E2=10,niu2=5.

Best Answer

There are two functions named diff().
One of them is part of the symbolic toolbox, and can do calculus differentiation of symbolic functions and symbolic expressions.
The other of them is regular MATLAB and does numeric differences on numeric arrays; this numeric one is more or less x(2:end) - x(1:end-1), just subtracting adjacent elements.
The form of call you used was more consistent with the symbolic version, but what you actually passed in was a function handle, then a particular numeric value for the current time, and then a numeric constant. This combination of parameters does not match either function.
I recommend that for generality you re-code as
function[dxdt]=sistemul_cu_ec_de_gradul_1(t,x)
   dxdt_2=x(1);
   dxdt_1=-(p1/p2)*x(1)-(p0/p2)*x(2)+(q2/p2)*d2_epsilon(t)+(q1/p2)*d_epsilon(t)+(q0/p2)*epsilon(t);
   %dxdt_1=-(p1/p2)*x(1)-(p0/p2)*x(2)+(q2/p2)*0+(q1/p2)*0+(q0/p2)*epsilon(t);
   dxdt=[dxdt_2; dxdt_1];
 end
function d2 = d2_epsilon(t)
  d2 = 0;
end
function d = d_epsilon(t)
  d = 0;
end
and manually update d2_epsilon and d_epsilon when you change epsilon.
If you have a need to not manually update the derivatives, then you will need to involve the symbolic toolbox before you call ode45, invoking
syms T
sym_epsilon = epsilon(T);
clear d_epsilon d2_epsilon
matlabFunction( diff(sym_epsilon, T), 'vars', T, 'file', 'd_epsilon.m');
matlabFunction( diff(sym_epsilon, T, 2), 'vars', T, 'file', 'd2_epsilon.m');
and not coding the all-zero functions that I showed.
Related Question