Problem 9.8 part A

I tried for number of time steps, for each time step found the largest eigenvalue for A using the power method. No Eigenvalue was found that is <= 1, so spectral radius, defined as magnitude of largest eigenvalue is not <=1 , hence system is unconditionally unstable.

One annoying thing with matlab, is that when I print a value it says it is 1, but when I find the difference between 1 and the value I do not always get a zero as expected. Even though I set format to long g. This is what the output below says that the largest eigenvalue is 1. using eps I see the largest eigenvalue is within the floating point accuracy from 1. so, is it 1 or not 1 ?

» help EPS

EPS Floating point relative accuracy.

EPS returns the distance from 1.0 to the next largest

floating point number. EPS is used as a default tolerance by

PINV and RANK, as well as several other MATLAB functions.

See also REALMAX, REALMIN.

» eps

2.22044604925031e-016

OUTPUT

» nma_problem_9_8_part_a

This solves problem 9.8 part a. finds if a matrix

is stable or not using the power method

Nasser Abbasi, April 29, 2002.

Enter N, the number of grid points (suggest 61):61

Enter n, number of iterations to find eigenvalues (suggest 10):10

h is 0.016949

trying for time steps from 0.000100 to 10.000000 at increment of 0.010000

largest eigenvalue after trying all time steps is 1

(maxLambda - 1) = 2.22045e-016

system is always unstable, no spectral radius found <= 1 for the time steps tried

Code

function nma_problem_9_8_part_a()

%

% This solves problem 9.8 part a. finds if a matrix

% is stable or not using the power method

%

% Nasser Abbasi, April 29, 2002.

clear all; help nma_problem_9_8_part_a;

N= input('Enter N, the number of grid points (suggest 61):');

%accuracy = input('Enter accuracy need to stop iteration (suggest 1e-8):');

%I first tried the loop below untill some user specified accuracy is reached

%for the v vector. This did not make any difference on the result. So I changed

%the program to force it to go to user specified number of multiplications

n = input('Enter n, number of iterations to find eigenvalues (suggest 10):');

L = 1;

c = 1;

h = L/(N-2);

fprintf('h is %f\n',h);

B = makeB(N);

I = eye(N);

x = ones(N,1);

firstTau = 0.0001;

lastTau = 10;

tauInc = 0.01;

fprintf('trying for time steps from %f to %f at increment of %f\n',firstTau,lastTau,tauInc);

maxLambda=-inf;

for(tau= firstTau: tauInc : lastTau)

A = (I - (c*tau/2*h) * B);

v = A^n*x;

v = v/norm(v);

mv= A*v;

lambda = mv(1)/v(1);

if(lambda > maxLambda )

maxLambda = lambda;

end

end

fprintf('largest eigenvalue after trying all time steps is %g\n',maxLambda);

fprintf('(maxLambda - 1) = %g\n',maxLambda-1);

if( maxLambda <= 1 )

fprintf('system can be stable for specific Tau, found a lambda less than one\n');

else

fprintf('system is always unstable, no spectral radius found <= 1 for the time steps tried\n');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function B=makeB(N)

B=zeros(N);

for(i=2:N-1)

for(j=1:N)

if( j == i+1)

B(i,j) = 1;

end

if( j == i-1)

B(i,j) = -1;

end

end

end

% now do the first and last rows

B(1,2) = 1;

B(1,N) = -1;

B(N,1) = 1;

B(N,N-1) = -1;