Ode45 and diff function problem

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.

Answer 1

[Walter Roberson]

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.

[Milu Mihai]

Ok i've modified what you said but now i got another error:

Undefined function 'd2_epsilon' for input arguments of type 'double'.

Error in serie_maxwell_kelvin/sistemul_cu_ec_de_gradul_1 (line 35)
   dxdt_1=-(p1/p2)*x(1)-(p0/p2)*x(2)+(q2/p2)*d2_epsilon(t)+(q1/p2)*d_epsilon(t)+(q0/p2)*epsilon(t);

Error in odearguments (line 88)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 114)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in serie_maxwell_kelvin (line 18)
[t,x]=ode45(@sistemul_cu_ec_de_gradul_1,t,[sigma_initial
dsigmadt_initial]);

The program modified it's the fallowing:

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');

   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');

   function [z] = epsilon(t)
  z=6;
  end

   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

   end

[Walter Roberson]

The symbolic stuff needs to be done before you call ode45.

By the way, you will want to change your initial conditions. With that epsilon function and those initial conditions, the result is all zeros.

[Milu Mihai]

Thanks alot,i noticed just now you actually said to call them before ode45.Made it work,thanks!