MATLAB: Index exceeds matrix dimensions.

matrix

clc; 
clear all; 
close all;
g = 1;
m = 1;
M = 1;
T = 1;
fun = @(x,y,V,T) ((2*pi*g*(1-(V.^2)).^2)/m^2).*y.*(((x-(m*V./sqrt(1-V.^2))).^2./(1-(V).^2))+y.^2).*exp(-(sqrt(m^2+x.^2+y.^2)-M)/T);
QQ_pl = @(V,T) integral2(@(y,x) fun(x,y,V,T), 0, inf, -inf, m*V./sqrt(1-V.^2));
QQ_mi = @(V,T) integral2(@(y,x) fun(x,y,V,T), 0, inf, m*V./sqrt(1-V.^2), inf);
V(1)=0.1
L=0;
dL=0.01;
t(1)=L+dL;
dV(1)=(QQ_mi(V(1),T(1))-QQ_pl(V(1),T(1)))*t(1)
for i=2:50
    t(i)=L+i*dL;
    dV(i)=(QQ_mi(V(i-1),T(i-1))-QQ_pl(V(i-1),T(i-1)))*t(i)
    V(i)=V(i-1)-dV(i)
end
plot(t,V)

Best Answer

You define T as a scalar

T = 1

but then in your for loop try to index out values from it like it is a vector

dV(i)=(QQ_mi(V(i-1),T(i-1))-QQ_pl(V(i-1),T(i-1)))*t(i)

You will therefore get this error message any time i>2 (since index is T(i-1)).

Just like you add values to V in each loop, you likely need to do something similar for T.

Related Solutions

MATLAB: Eig() gives wrong results for complex eigenvalues in comparison to eigs()

Hi Lennart,

It seems there are a few things that are not correct with the last part of the code:

1. Instead of li(qq), it should be li(qq,qq).

2. You should not use (.'), instead (') should be used. (.') computes the nonconjugate transpose of a matrix whereas (.) computes complex conjugate transpose. Since the eigenvalues and eigenvectors could be complex in nature even for a real matrix, the latter is the correct form to be used.

3. The code below gives you the desired results:

for qq = 1:length(L1)
  %





  right1(:, qq) = (A - L1(qq,qq) * B) * X1(:, qq);
  left1(qq, :) = Y1(:, qq)' * (A - L1(qq,qq) * B);
  %
  right2(:, qq) = (A - L21(qq,qq) * B) * X2(:, qq);
  left2(qq, :) = Y2(:, qq)' * (A - L22(qq,qq)' * B);
  %
  right3(:, qq) = (A - L3(qq,qq) * B) * X3(:, qq);
  left3(qq, :) = Y3(:, qq)' * (A - L3(qq,qq) * B);
  %
  right4(:, qq) = (A - L41(qq,qq) * B) * X4(:, qq);
  left4(qq, :) = Y4(:, qq)' * (A - L42(qq,qq)' * B);
  %
  right5(:, qq) = (A - L5(qq,qq) * B) * X5(:, qq);
  left5(qq, :) = Y5(:, qq)' * (A - L5(qq,qq) * B);
  %
  right6(:, qq) = (A - L61(qq,qq) * B) * X6(:, qq);
  left6(qq, :) = Y6(:, qq)' * (A - L62(qq,qq)' * B);
end
[left1 right1;
left2 right2;
left3 right3;
left4 right4;
left5 right5;
left6 right6]

Please note that I have changed Li2(qq,qq) to Li2(qq,qq)' for the computation of left2, left4 and left6 in the 2nd, 4th and 6th case. Here is why: As per the documentation of eig,

[Y2, L22] = eigs(A', B');

ensures that the solution satisfies the following:

A'*Y2 = B'*Y2*L22   % which gives  
Y2'*A = L22'*Y2'*B  % (taking transpose of both sides)
Y2'*A - L22'*Y2'*B = 0  
Y2'*(A-L22'*B) = 0

(since values of L22 are taken one by one in your calculations, so it can behave as a scalar and moving its position does not change the result). However, taking its (complex conjugate) transpose is crucial since it is a complex number. Instead of computing the residuals inside a loop, it can be written as one single statement (treating L22 as a matrix):

left2 = Y2'*A - L22'*Y2'*B;

The same reasoning applies for cases 4 and 6. Please note this reasoning does not apply to cases 1, 3 and 5 and the original form of equation stays valid.

I hope the above helps.

Nalini

MATLAB: How to display velocity for each loop at 20 seconds

You have an error in your code:

fprint('%.0m/s\n', vj(,2)

This code should work for you:

Time=20; %Time the jumper free-falls
Cd=[0.25,.50,.75]; %varying coefficients of drag in kg/m
m=[40,60,80]; %varying masses in kg
g=9.81; %accleration of jumper due to gravity
n=100; %number steps
deltaT=Time/(n-1);
v(1)=0; %initialVelocity
t=linspace(0,Time,n); %the time interval of the jumper
for j=1:3
  hold on
  for i=2:n
    v(i)=v(i-1)+(g-(Cd(j)/m(2))*v(i-1)*abs(v(i-1)))*deltaT;
  end 
  vj=v;
  plot(t,vj)
end
fprintf('%.0f m/s\n',vj(1))
fprintf('%.0f m/s\n',vj(2))

Best Answer

Related Solutions

MATLAB: Eig() gives wrong results for complex eigenvalues in comparison to eigs()

MATLAB: How to display velocity for each loop at 20 seconds

Related Question