It is possible to update a matrix each step?
7 views (last 30 days)
Show older comments
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 Comments
Babak
on 7 May 2013
Yes, the variable K is updated according to K1, K2, K3 every step in the for loop with index i.
Accepted Answer
James Kristoff
on 7 May 2013
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 Comments
More Answers (0)
See Also
Categories
Find more on Matrix Indexing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!