Talk:Forum for GNU Octave: Difference between revisions

Jump to navigation Jump to search
Delete question give hint where to find help.
(why output always zero after each time integration points?)
 
(Delete question give hint where to find help.)
Tag: Replaced
Line 1: Line 1:
Dear All,
Please ask questions at Octave's mailing-list, see [[FAQ#Where_do_I_get_additional_help.3F]].
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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
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

Navigation menu