https://wiki.octave.org/wiki/api.php?action=feedcontributions&user=95.65.232.248&feedformat=atomOctave - User contributions [en]2024-03-29T05:20:26ZUser contributionsMediaWiki 1.39.2https://wiki.octave.org/wiki/index.php?title=Talk:Forum_for_GNU_Octave&diff=12851Talk:Forum for GNU Octave2020-04-14T13:09:56Z<p>95.65.232.248: why output always zero after each time integration points?</p>
<hr />
<div>Dear All,<br />
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?<br />
for the main .m file it is as below:<br />
<br />
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br />
alpha = [0.0020*2*pi 8*pi 0];<br />
tspan = [0 4] ;<br />
frunb=zeros(size(K,1),1); %%% unbalance vector<br />
frunb(5,1) = 1; %%% <br />
nr = 3;<br />
options = odeset;<br />
[t,y] = ode45(@deriv,tspan,zeros(2*nr,1),options,M,K,C,G,alpha,frunb);<br />
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br />
<br />
<br />
and below is the function deriv (saved in a different file named as deriv.m)<br />
<br />
<br />
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br />
function dz = deriv(t,z,M,K,C,G,alpha,ubforce)<br />
<br />
phi = alpha(1)*t*t + alpha(2)*t + alpha(3);<br />
d_phi = 2*alpha(1)*t + alpha(2);<br />
dd_phi = 2*alpha(1);<br />
A = [(-inv(M)*(C-j*phi*G)) (-inv(M)*K); eye(size(M)) zeros(size(M))];<br />
[U,D] = eig(A);<br />
uncoupled_A = (inv(U))*A*(U);<br />
B = [inv(M)*ubforce;zeros(size(M))*ubforce];<br />
uncoupled_B = inv(U) * B;<br />
double_uncoupled_A_1 = [(uncoupled_A(1,1) + uncoupled_A(2,2)) (-uncoupled_A(1,1)*uncoupled_A(2,2)); 1 0]; <br />
double_uncoupled_A_2 = [(uncoupled_A(3,3) + uncoupled_A(4,4)) (-uncoupled_A(3,3)*uncoupled_A(4,4)); 1 0];<br />
double_uncoupled_A_3 = [(uncoupled_A(5,5) + uncoupled_A(6,6)) (-uncoupled_A(5,5)*uncoupled_A(6,6)); 1 0];<br />
U11 = [uncoupled_A(1,1) uncoupled_A(2,2); 1 1];<br />
U22 = [uncoupled_A(3,3) uncoupled_A(4,4); 1 1];<br />
U33 = [uncoupled_A(5,5) uncoupled_A(6,6); 1 1];<br />
U = U(1:6,1:6);<br />
U1 = blkdiag(U11,U22,U33);<br />
double_uncoupled_A = blkdiag(double_uncoupled_A_1,double_uncoupled_A_2,double_uncoupled_A_3);<br />
double_uncoupled_B = U1 * inv(U) * B(1:6,:);<br />
gain = [1 0 0 0 0 0] * (phi^2)*exp(j*phi*t)*z;<br />
dz = (double_uncoupled_A * z) + double_uncoupled_B * gain';<br />
endfunction</div>95.65.232.248