It is possible to update a matrix each step?

조회 수: 7 (최근 30일)
Maikel
Maikel 2013년 5월 7일
Hi everyone, I'm trying to program a matrix dependant of the values of certain variables that depend on time (matrix K in the code below). I've programmed but I really don't know if the program is taking the first value that matches the condition (it compares x values) or it evaluates each step the values of K1, K2, K3.
I've tried before with an K(:,i+1) but it seems that it doesn't work with matrices (or I don't know how)
Please, would you give me some advice?
Thx in advance. PS: I'm Spanish please excuse me for my poor english
function [K] = Newm3
clear;clc;
h=1;
x=[0;0;0];
v=[sqrt(2*9.81*h);0;0];
F=[0;0;0];
M=[600 0 0;0 174.875 0;0 0 321.186];
C=[0 0 0;0 0 0;0 0 0];
K=[1e7 -1e7 0;-1e7 1e7+8.022e7 -8.022e7;0 -8.022e7 8.024e7];
a=(M^-1)*(F-C*v-K*x);
al=1/6;
be=1/2;
% Time step
ti = 0. ;
tf = 1. ;
dt = 0.0005 ;
t = ti:dt:tf ;
nt = fix((tf-ti)/dt) ;
n = length(M) ;
for i = 1:nt
part1=M*((1/(al*dt^2))*x(:,i)+(1/(al*dt))*v(:,i)+((1/(2*al))-1)*a(:,i));
part2=C*((be/(al*dt))*x(:,i)+((be/al)-1)*v(:,i)+((be/al)-2)*(dt/2)*a(:,i));
x(:,i+1)=((1/(al*(dt)^2))*M+(be/(al*dt))*C+K)^-1*...
(F+part1+part2);
if x(1,i+1)-x(2,i+1) > 0,
K1=1e8;
else
K1=0;
end
dx=abs(x(2,i+1)-x(3,i+1));
if dx<=3.848e-5,
K2=8.018e7;
elseif 3.848e-5<dx<1.103e-4
K2=-1.295e8;
elseif 1.103e-4<dx<1.821e-4
K2=4.618e4;
else
K2=342.036;
end
if x(3,i+1)-x(3,i)<=0
K3=1.766e4;
elseif x(3,i+1)-x(3,i)>0
K3=3.208e3;
else
K3=56.104
end
K=[K1 -K1 0;-K1 K1+K2 -K2;0 -K2 K2+K3];
a(:,i+1)=(1/(al*dt^2))*(x(:,i+1)-x(:,i))-(1/(al*dt))*v(:,i)-((1/(2*al))-1)*a(:,i);
v(:,i+1)=v(:,i)+(1-be)*dt*a(:,i)+be*dt*a(:,i+1);
end
figure;
hold on
d1=plot(t,x(1,:)) ;
set(d1,'Color','red')
d2=plot(t,x(2,:)) ;
set(d2,'Color','green')
d3=plot(t,x(3,:)) ;
set(d3,'Color','blue')
hold off
end
  댓글 수: 6
Babak
Babak 2013년 5월 7일
Yes, the variable K is updated according to K1, K2, K3 every step in the for loop with index i.
Maikel
Maikel 2013년 5월 8일
Thank you very much, I really appreciate your comments.

댓글을 달려면 로그인하십시오.

채택된 답변

James Kristoff
James Kristoff 2013년 5월 7일
It is possible to update a matrix in each iteration of a loop. Your code is currently updating the matrix K every iteration and then using the updated values in the next loop. If you want to capture every version of this matrix, you could do something similar to the following:
myMatrix = [0,0; 0,0];
for i = 1:10
myMatrix(:,:,i) = [i, i+1; i^2, i/2];
end
The : operator means, all values in this dimension, and since you want to create a copy of the matrix, which already has 2 dimensions, so you need to make the iterator the third dimension.
I made some comments about the code you posted below:
function [K] = Newm3
% CLEAR is not needed inside a function like this. Functions have their
% own workspace, which initially contains only the inputs to the function,
% unless the function contains persistent variables.
%
% clear;clc;
%
clc;
% try to name your variables something that makes sense
h=1; %height?
x=[0;0;0]; % position?
v=[sqrt(2*9.81*h);0;0]; % velocity? if so, I believe this equation is incorrect.
F=[0;0;0]; % force?
M=[600 0 0;...
0 174.875 0;...
0 0 321.186]; % Mass?
% this can also be done with C = zeros(3);
C=[0 0 0;...
0 0 0;...
0 0 0]; % dampener coeffiecent matrix?
K=[ 1e7 -1e7 0;...
-1e7 1e7+8.022e7 -8.022e7;...
0 -8.022e7 8.024e7]; % stiffness coeffiecient matrix?
a=(M^-1)*(F-C*v-K*x); % accelleration
al=1/6; %?
be=1/2; %?
% Time step
ti = 0.;
tf = 1.;
dt = 0.0005;
t = ti:dt:tf;
nt = fix((tf-ti)/dt);
% unused variable
% n = length(M);
for i = 1:nt
part1 = M*((1/(al*dt^2))*x(:,i) + (1/(al*dt))*v(:,i) + ((1/(2*al))-1)*a(:,i));
part2 = C*((be/(al*dt)) *x(:,i) + ((be/al)-1)*v(:,i) + ((be/al)-2)*(dt/2)*a(:,i));
% F starts as [0;0;0] and never changes.
x(:,i+1) = ((1/(al*(dt)^2))*M + (be/(al*dt))*C+K)^-1*(F+part1+part2);
if x(1,i+1)-x(2,i+1) > 0,
K1=1e8;
else
K1=0;
end
dx=abs(x(2,i+1)-x(3,i+1));
if dx<=3.848e-5,
K2=8.018e7;
elseif (3.848e-5<dx) && (dx<1.103e-4) % the previous syntax is invalid
K2=-1.295e8;
elseif (1.103e-4<dx) && (dx<1.821e-4)
K2=4.618e4;
else
K2=342.036;
end
if x(3,i+1)-x(3,i)<=0
K3=1.766e4;
elseif x(3,i+1)-x(3,i)>0
K3=3.208e3;
else
K3=56.104;
end
K=[ K1 -K1 0;...
-K1 K1+K2 -K2;...
0 -K2 K2+K3];
a(:,i+1) = (1/(al*dt^2))*(x(:,i+1) - x(:,i))-(1/(al*dt))*v(:,i) - ((1/(2*al))-1)*a(:,i);
v(:,i+1) = v(:,i) + (1-be)*dt*a(:,i) + be*dt*a(:,i+1);
end
figure;
hold on
d1=plot(t,x(1,:)) ;
set(d1,'Color','red')
d2=plot(t,x(2,:)) ;
set(d2,'Color','green')
d3=plot(t,x(3,:)) ;
set(d3,'Color','blue')
hold off
end

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Matrix Indexing에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by