Talk:Forum for GNU Octave: Difference between revisions

From Octave
Jump to navigation Jump to search
(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

Revision as of 03:08, 15 April 2020

Please ask questions at Octave's mailing-list, see FAQ#Where_do_I_get_additional_help.3F.

Retrieved from "https://wiki.octave.org/wiki/index.php?title=Talk:Forum_for_GNU_Octave&oldid=12852"