Talk:Forum for GNU Octave

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

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