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 个评论
显示 4更早的评论
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 Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by