Thread Subject:

Command window output

I'm in a numerical analysis class and we've been going over Jacobi iteration and Gauss-Seidel. I understand the problems, and how to code them, but i'm having a problem with unwanted answers being shown in the command window. As i'm using the diary command to print the code and the answers, this is a problem. Here's my code:

%% P2.78
a1 = [7 -3 4];
a2 = [2 5 3];
a3 = [-3 2 6];
b = [6 -5 2];

x1(1)=0;
x2(1)=0;
x3(1)=0;
i(1)=1;

epsilon=1.*10^-5;
residual=10;


while residual > epsilon
    x1(i+1)=x1(i) + (1./a1(1)).*(b(1) - a1(1).*x1(i) - a1(2).*x2(i) - a1(3).*x3(i));
    x2(i+1)=x2(i) + (1./a2(2)).*(b(2) - a2(1).*x1(i) - a2(2).*x2(i) - a2(3).*x3(i));
    x3(i+1)=x3(i) + (1./a3(3)).*(b(3) - a3(1).*x1(i) - a3(2).*x2(i) - a3(3).*x3(i));
    i=i+1;
    residual(i)=abs(7.*x1(i) - 3.*x2(i) +4.*x3(i) -6);
end

fid=fopen('P278.txt','w+');
fprintf(fid,'i x1 x2 x3 Residual\n');
for k=1:max(i);
    fprintf(fid,'%u %f %f %f %f\n',k,x1(k),x2(k),x3(k),residual(k));
end
fclose(fid);
type('P278.txt')


%% P2.81
% Part a
a1 = [4.63 -1.21 3.22];
a2 = [-3.07 5.48 2.11];
a3 = [1.26 3.11 4.57];
b = [2.22 -3.17 5.11];

x1(1)=0;
x2(1)=0;
x3(1)=0;
i=1;

epsilon=1.*10^-5;
residual=10;

while residual > epsilon
    x1(i+1)=x1(i) + (1./a1(1)).*(b(1) - a1(1).*x1(i) - a1(2).*x2(i) - a1(3).*x3(i));
    x2(i+1)=x2(i) + (1./a2(2)).*(b(2) - a2(1).*x1(i) - a2(2).*x2(i) - a2(3).*x3(i));
    x3(i+1)=x3(i) + (1./a3(3)).*(b(3) - a3(1).*x1(i) - a3(2).*x2(i) - a3(3).*x3(i));
    i=i+1;
    residual(i)=abs(4.63.*x1(i) - 1.21.*x2(i) +3.22.*x3(i) -2.22);
end

fid=fopen('P281a.txt','w+');
fprintf(fid,'i x1 x2 x3 Residual\n');
for k=1:max(i);
    fprintf(fid,'%u %f %f %f %f\n',k,x1(k),x2(k),x3(k),residual(k));
end
fclose(fid);
type('P281a.txt')

%Part b
a1 = [4.63 -1.21 3.22];
a2 = [-3.07 5.48 2.11];
a3 = [1.26 3.11 4.57];
b = [2.22 -3.17 5.11];

x1(1)=0;
x2(1)=0;
x3(1)=0;
i=1;

epsilon=1.*10^-5;
residual=10;

while residual > epsilon
    x1(i+1)=x1(i) + (1./a1(1)).*(b(1) - a1(1).*x1(i) - a1(2).*x2(i) - a1(3).*x3(i));
    x2(i+1)=x2(i) + (1./a2(2)).*(b(2) - a2(1).*x1(i+1) - a2(2).*x2(i) - a2(3).*x3(i));
    x3(i+1)=x3(i) + (1./a3(3)).*(b(3) - a3(1).*x1(i+1) - a3(2).*x2(i+1) - a3(3).*x3(i));
    i=i+1;
    residual(i)=abs(4.63.*x1(i) - 1.21.*x2(i) +3.22.*x3(i) -2.22);
end

fid=fopen('P281b.txt','w+');
fprintf(fid,'i x1 x2 x3 Residual\n');
for k=1:max(i);
    fprintf(fid,'%u %0.5f %0.5f %0.5f %0.5f\n',k,x1(k),x2(k),x3(k),residual(k));
end
fclose(fid);
type('P281b.txt')

When i run the code, i get the tables i was trying to make with fprintf, but i also get ans=63,ans=59,etc. Where is this coming from?

Nevermind, i just had to restart matlab and it took care of it.

Jeremy wrote:

 > When i run the code, i get the tables i was trying to make with fprintf,
 > but i also get ans=63,ans=59,etc. Where is this coming from?

I do not happen to notice why that is happening, but I happened to notice an
inefficiency in your code:


> epsilon=1.*10^-5;
> residual=10;
>
>
> while residual > epsilon
> x1(i+1)=x1(i) + (1./a1(1)).*(b(1) - a1(1).*x1(i) - a1(2).*x2(i) -
> a1(3).*x3(i));
> x2(i+1)=x2(i) + (1./a2(2)).*(b(2) - a2(1).*x1(i) - a2(2).*x2(i) -
> a2(3).*x3(i));
> x3(i+1)=x3(i) + (1./a3(3)).*(b(3) - a3(1).*x1(i) - a3(2).*x2(i) -
> a3(3).*x3(i));
> i=i+1;
> residual(i)=abs(7.*x1(i) - 3.*x2(i) +4.*x3(i) -6);
> end

In this loop and the other loops like it, you grow the array 'residual' on
each trip through the loop. It would be more efficient if you could
pre-allocate residual, even if only as an (over-?) estimate of the number of
steps you will take. You might still end up growing the array if you
underestimate, but the growth will be postponed until the size you pre-allocated.

As you grow the array, the while test you have,

while residual > epsilon

becomes equivalent, for Matlab purposes, to

while all(residual > epsilon)

As you eventually get to a small enough residual somewhere in the array, that
will happen to work, but it is inefficient. If you stick with the structure
where you grow residual on every trip through, then you could increase the
efficiency by coding

while residual(end) > epsilon

Alternately, another version that would work even if you pre-allocate, would be

while residual(i) > epsilon


Also note that each trip through the loops, you re-calculate the
subexpressions (1./a1(1)) and (1./a2(2)) and (1./a3(3)) . That suggests that
you should move that outside the loops:

oneovera1 = 1./a1(1);
oneovera2 = 1./a2(2);
oneovera3 = 1./a3(3);

then you would have lines such as:

x1(i+1)=x1(i) + ...
oneovera1.*(b(1) - a1(1).*x1(i) - a1(2).*x2(i) - a1(3).*x3(i));

I don't know how many trips through the loop are typical, but removing three
divisions per cycle would surely be an improvement.

Walter Roberson <roberson@hushmail.com> wrote in message <hm703a$f5i$1@canopus.cc.umanitoba.ca>...
> Jeremy wrote:
>
> > When i run the code, i get the tables i was trying to make with fprintf,
> > but i also get ans=63,ans=59,etc. Where is this coming from?
>
> I do not happen to notice why that is happening, but I happened to notice an
> inefficiency in your code:
>
>
> > epsilon=1.*10^-5;
> > residual=10;
> >
> >
> > while residual > epsilon
> > x1(i+1)=x1(i) + (1./a1(1)).*(b(1) - a1(1).*x1(i) - a1(2).*x2(i) -
> > a1(3).*x3(i));
> > x2(i+1)=x2(i) + (1./a2(2)).*(b(2) - a2(1).*x1(i) - a2(2).*x2(i) -
> > a2(3).*x3(i));
> > x3(i+1)=x3(i) + (1./a3(3)).*(b(3) - a3(1).*x1(i) - a3(2).*x2(i) -
> > a3(3).*x3(i));
> > i=i+1;
> > residual(i)=abs(7.*x1(i) - 3.*x2(i) +4.*x3(i) -6);
> > end
>
> In this loop and the other loops like it, you grow the array 'residual' on
> each trip through the loop. It would be more efficient if you could
> pre-allocate residual, even if only as an (over-?) estimate of the number of
> steps you will take. You might still end up growing the array if you
> underestimate, but the growth will be postponed until the size you pre-allocated.
>
> As you grow the array, the while test you have,
>
> while residual > epsilon
>
> becomes equivalent, for Matlab purposes, to
>
> while all(residual > epsilon)
>
> As you eventually get to a small enough residual somewhere in the array, that
> will happen to work, but it is inefficient. If you stick with the structure
> where you grow residual on every trip through, then you could increase the
> efficiency by coding
>
> while residual(end) > epsilon
>
> Alternately, another version that would work even if you pre-allocate, would be
>
> while residual(i) > epsilon
>
>
> Also note that each trip through the loops, you re-calculate the
> subexpressions (1./a1(1)) and (1./a2(2)) and (1./a3(3)) . That suggests that
> you should move that outside the loops:
>
> oneovera1 = 1./a1(1);
> oneovera2 = 1./a2(2);
> oneovera3 = 1./a3(3);
>
> then you would have lines such as:
>
> x1(i+1)=x1(i) + ...
> oneovera1.*(b(1) - a1(1).*x1(i) - a1(2).*x2(i) - a1(3).*x3(i));
>
> I don't know how many trips through the loop are typical, but removing three
> divisions per cycle would surely be an improvement.

The loops typically only run through ~250 iterations at max. If you start using over-relaxation or under-relaxation and use bad values for w, then it could get up into the thousands, which still wouldn't make a noticeable difference in computing time for my purposes.

Thank you for the suggestions though. I really need to start looking at how to improve the efficiency of my coding, and now's a pretty good time to start i guess. Would you suggest i pre-allocate my x1, x2, and x3 array also?