why output always zero after each time integration points?

Dear All,

I am trying to write an transient analysis solver for a multi degree of freedom system. it consists of mass, stiffness, damping, gyroscopic damping matrices and so on. my initial condition is declared by zeros(2*nr,1). after I run the code, the output y which is a matrix of 6x41 has all it's rows and columns just zero. It would be really apprectiated if some one has an advice / opinion about what am I doing wrong?

for the main .m file it is as below:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alpha = [0.0020*2*pi 8*pi 0];
tspan = [0 4] ;
frunb=zeros(size(K,1),1); %%% unbalance vector
frunb(5,1) = 1; %%% 
nr = 3;
options = odeset;
[t,y] = ode45(@deriv,tspan,zeros(2*nr,1),options,M,K,C,G,alpha,frunb);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

and below is the function deriv (saved in a different file named as deriv.m)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dz = deriv(t,z,M,K,C,G,alpha,ubforce)
  


  phi = alpha(1)*t*t + alpha(2)*t + alpha(3);
  d_phi = 2*alpha(1)*t + alpha(2);
  dd_phi = 2*alpha(1);

  A = [(-inv(M)*(C-j*phi*G)) (-inv(M)*K); eye(size(M)) zeros(size(M))];
  [U,D] = eig(A);
  uncoupled_A = (inv(U))*A*(U);

  B = [inv(M)*ubforce;zeros(size(M))*ubforce];
  uncoupled_B = inv(U) * B;

  double_uncoupled_A_1 = [(uncoupled_A(1,1) + uncoupled_A(2,2)) (-uncoupled_A(1,1)*uncoupled_A(2,2)); 1 0]; 
  double_uncoupled_A_2 = [(uncoupled_A(3,3) + uncoupled_A(4,4)) (-uncoupled_A(3,3)*uncoupled_A(4,4)); 1 0];
  double_uncoupled_A_3 = [(uncoupled_A(5,5) + uncoupled_A(6,6)) (-uncoupled_A(5,5)*uncoupled_A(6,6)); 1 0];

  U11 = [uncoupled_A(1,1)  uncoupled_A(2,2); 1 1];
  U22 = [uncoupled_A(3,3)  uncoupled_A(4,4); 1 1];
  U33 = [uncoupled_A(5,5)  uncoupled_A(6,6); 1 1];

  U = U(1:6,1:6);
  U1 = blkdiag(U11,U22,U33);
  double_uncoupled_A = blkdiag(double_uncoupled_A_1,double_uncoupled_A_2,double_uncoupled_A_3);

  double_uncoupled_B = U1 * inv(U) * B(1:6,:);

  gain = [1 0 0 0 0 0] * (phi^2)*exp(j*phi*t)*z;
  dz = (double_uncoupled_A * z) + double_uncoupled_B * gain';


endfunction

[Walter Roberson]

You should be asking in an Octave resource for Octave code.

[kadir]

I wrote this code in my own personal pc. however I run it on my pc at the office. due to home office situation right now, I am working at home and no access to my pc at office.

[Walter Roberson]

Yes, but we here are not qualified to talk about potential problems in Octave's ode45 code, or what oddities might exist in using inv() with it, or so on.

By the way, never use inv() for these kinds of purposes. inv(M)*ubforce is better written as M\ubforce and likewise -inv(M)*K is better written as M\K .

If you think it is important to use inv() for your situation, then at least assign it to a variable so that you do not keep recomputing it.

And see

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

The details of passing extra parameters for Octave's ode45 are probably different than the details of passing extra paramters to MATLAB's ode45(). The ability to pass extra parameters that way has not been documented for a good 15 years for MATLAB.