im having a lot of trouble with my project code, and i dont know very much about matlab besides the codes he's given us in this class and what ive been trying to figure out all day today, i think my errors are probably simple for someone who is more familar with, im just somehow not defining my 3d array correctly.

can anyone look at it and help

12/13/2008 6:59:57 PM

what specifically is the problem? what are your error messages? i'm not too familiar with 3d arrays, but i have a pretty good handle on general matlab syntax.

12/13/2008 7:22:27 PM

pcg3d is my main program in it , it calls pcgssor- both seem to be having trouble witht he same arrary "up"

I can email u the code, pm me an email address
thank u for any help, my professor is out of town or something

>> pcg3d
??? Attempted to access up(2,2,2); index out of bounds because
size(up)=[21,21,1].

Error in ==> pcg3d at 28
ar(i,j,l) = -(cond(up(i,j,l))+cond(up(i,j,l+1)))*.5;

>> edit cond.m
>> pcgssor
??? Input argument "up" is undefined.

Error in ==> pcgssor at 5
u = up;

12/13/2008 7:28:20 PM

Your 3rd dimension of the array "up" looks like it only has 1 value, but your statement at 28 is trying to step up to a second value, that doesn't exist.

Either you're not initializing your 3rd dimension correctly, or your algorithm to generate your "ar" array is wrong.

12/13/2008 7:42:24 PM

>> pcg3d
??? Attempted to access up(2,2,2); index out of bounds because
size(up)=[21,21,1].

12/13/2008 7:43:20 PM

clear;
% This progran solves -(K(u)ux)x - (K(u)uy)y – (K(u)uz)z = f.
% K(u) is defined in the function cond(u).
% The Picard nonlinear method is used.
% The solve step is done in the subroutine cgssor.
% It uses the PCG method with SSOR preconditioner.
maxmpic = 50;
tol = .001;
n = 20;
up = zeros(n+1);
rhs = zeros(n+1);
up = zeros(n+1);
h = 1./n;
% Defines the right side of PDE.
for l=2:n
for j = 2:n
for i = 2:n
rhs(i,j,l) = h*h*h*200.*sin(3.14*(i-1)*h)*sin(3.14*(j-1)*h)*sin(3.14*(l-1)*h);
end
end
end
% Start the Picard iteration.
for mpic=1:maxmpic
% Defines the 7 nonzero components in the 3d array?.
for l = 2:n
for j = 2:n
for i = 2:n
ar(i,j,l) = -(cond(up(i,j,l))+cond(up(i,j,l+1)))*.5;
ag(i,j,l) = -(cond(up(i,j,l))+cond(up(I,j,l-1)))*.5;
an(i,j,l) = -(cond(up(i,j,l))+cond(up(i,j+1,l)))*.5;
as(i,j,l) = -(cond(up(i,j,l))+cond(up(i,j-1,l)))*.5;
ae(i,j,l) = -(cond(up(i,j,l))+cond(up(i+1,j,l)))*.5;
aw(i,j,l) = -(cond(up(i,j,l))+cond(up(i-1,j,l)))*.5;
ac(i,j,l) = -(ar(i,j,l)+ ag(i,j,l)+ an(i,j,l)+as(i,j,l)+ ae(i,j,l)+aw(i,j,l));
end
end
end
%
% The solve step is done by PCG with SSOR preconditioner.
%
[u , mpcg] = pcgssor(ar,ag,an,as,aw,ae,ac,up,rhs,n);
%
errpic = max(max(abs(up(2:n,2:n,2:n)-u(2:n,2:n,2:n))));
fprintf('Picard iteration = %6.0f\n',mpic)
fprintf('Number of PCG iterations = %6.0f\n',mpcg)
fprintf('Picard error = %6.4e\n',errpic)
fprintf('Max u = %6.4f\n', max(max(u)))
up = u;
if (errpic<tol)
break;
end
end
slice(u, 10,10,10) % creates color coded 3D plot
colorbar

12/13/2008 8:05:26 PM

it calls this one:
% PCG subroutine with SSOR preconditioner

function [u , mpcg]= pcgssor(ar,ag,an,as,aw,ae,ac,up,rhs,n)
w = 1.5;
u = up;
r = zeros(n+1);
rhat = zeros(n+1);
q = zeros(n+1);
p = zeros(n+1);
% Use the previous Picard iterate as an initial guess for PCG.
for l =2:n
for j = 2:n
for i = 2:n
r(i,j,l) = rhs(i,j,l)-(ac(i,j,l)*up(i,j,l) ...
+aw(i,j,l)*up(i-1,j,l)+ae(i,j,l)*up(i+1,j,l) ...
+as(i,j,l)*up(i,j-1,l)+an(i,j,l)*up(i,j+1,l)) ...
+ag(i,j,l)*up(i,j,l-1)+ar(i,j,l)*up(i,j,l+1);
end
end
end
error = 1. ;
m = 0;
rho = 0.0;
while ((error>.0001)&(m<200))
m = m+1;
oldrho = rho;
% Execute SSOR preconditioner.
for l= 2:n
for j= 2:n
for i = 2:n
rhat(i,j,l) = w*(r(i,j,l)-aw(i,j,l)*rhat(i-1,j,l) ...
-as(i,j,l)*rhat(i,j-1,l) ...
-ag(i,j,l)*rhat(i,j,l-1))/ac(i,j,l);
end
end
end
for l=2:n
for j= 2:n
for i = 2:n
rhat(i,j,l) = ((2.-w)/w)*ac(i,j,l)*rhat(i,j,l);
end
end
end
for l=-1:2
for j= n:-1:2
for i = n:-1:2
rhat(i,j,l) = w*(rhat(i,j,l)-ae(i,j,l)*rhat(i+1,j,l) ...
-an(i,j,l)*rhat(i,j+1,l))...
-ar(i,j,l)*rhat(i,j,l+1)/ac(i,j,l);
end
end
end
% Find conjugate direction.
rho = sum(sum(r(2:n,2:n,2:n).*rhat(2:n,2:n,2:n)));
if (m==1)
p = rhat;
else
p = rhat + (rho/oldrho)*p ;
end
% Execute 3d array product q = Ap.
for l =2:n
for j = 2:n
for i = 2:n
q(i,j,l)=ac(i,j,l)*p(i,j,l)+aw(i,j,l)*p(i-1,j,l) ...
+ae(i,j,l)*p(i+1,j,l)+as(i,j,l)*p(i,j-1,l) ...
+an(i,j,l)*p(i,j+1,l) +ag(i,j,l)*p(i,j,l-1)...
+ar(i,j,l)*p(I,j,l+1);
end
end
end
% Find steepest descent.
alpha = rho/sum(sum(p.*q));
u = u + alpha*p;
r = r - alpha*q;
error = max(max(abs(r(2:n,2:n,2:n))));
end
mpcg = m;

12/13/2008 8:06:43 PM

i think ar is being intialized fine but it does say "might be grown using indexing consider preallocating for speed"

perhaps i am not intialzing the thrid dimension correctly in the first few lines with up?
can either of u see whats wrong

this isnt a class in matlab and we didnt have to have a class in it as prereqs were just supposed to modify code he gives us so thats why im having trouble, hes usually very helpful but isnt responding.

12/13/2008 8:09:28 PM

up = zeros(n+1);
this statement defines a 2d array, and was causing problems. i changed it to up = zeros(21,21,21). that solution may not be practical for your assignment, but it worked.

also matlab is case sensitive, so you'll need to fix this:

ag(i,j,l) = -(cond(up(i,j,l))+cond(up(I,j,l-1)))*.5;


it craps out when it calls pcgssor

oh snap, more code. didn't see that. wait for it. . .

[Edited on December 13, 2008 at 8:23 PM. Reason : .]

12/13/2008 8:22:55 PM

two issues with your pcgssor function:

rhat=zeros(n+1)
again, this creates a 2d array, pcssor is looking for the third dimension. fixed it by making rhat=zeros(21,21,21);

also, further down you have this:
for l=-1:2
for j= n:-1:2
for i = n:-1:2

array indexes must be positive integers. i'm not sure what you intent here is.

oh, i think i see, are you trying to increment backwards down from n to 2?

[Edited on December 13, 2008 at 8:42 PM. Reason : .]

12/13/2008 8:34:23 PM

that was super helpful! i think i have just one more issue with moving from 2-d to 3-d
i guess row, p, q are not set up with 3d's correct or is it that i need to find another way to "divide these 3d arrays"

??? Error using ==> mrdivide
Input arguments must be 2-D.

Error in ==> pcgssor at 72
alpha = rho/sum(sum(p.*q));

12/13/2008 8:57:31 PM

yeah, there's something wierd going on there, i'm looking at it now. not sure what the deal is.

ok,
rho is (1x1x19)
p is (21x21x21)
q is (21x21x20)

the elements all need to be of the same number of dimensions. this could be accomplished explicitly by doing something like, but i don't know if it will give you answer you're looking for:
alpha = rho./sum(sum(p(1,1,1:19).*q(1,1,1:19)));

as an aside, a large part of these intermediate arrays are NaNs.

[Edited on December 13, 2008 at 9:11 PM. Reason : .]

12/13/2008 9:06:01 PM

without looking at any other code, i'm going to guess you need to do a component-wise divide (put a period in front of the division symbol).

12/13/2008 9:08:05 PM

this thread is proof i need a girlfriend, but anyway...

what exactly should the array up contain? as it stands, you are performing the cond on an array of zeros. this is what i'm referring to:

ar(i,j,l) = -(cond(up(i,j,l))+cond(up(i,j,l+1)))*.5;
ag(i,j,l) = -(cond(up(i,j,l))+cond(up(i,j,l-1)))*.5;
an(i,j,l) = -(cond(up(i,j,l))+cond(up(i,j+1,l)))*.5;
as(i,j,l) = -(cond(up(i,j,l))+cond(up(i,j-1,l)))*.5;
ae(i,j,l) = -(cond(up(i,j,l))+cond(up(i+1,j,l)))*.5;
aw(i,j,l) = -(cond(up(i,j,l))+cond(up(i-1,j,l)))*.5;

12/13/2008 9:50:56 PM

^^^ last I checked, component-wise operations don't care about matching dimensions, they just don't compute for mismatched components.

12/13/2008 9:55:12 PM

^wierd, i gives me an error if the dims don't match. i'm using R2006b, maybe they've updated it since then.

12/13/2008 10:01:38 PM

the period in front of the divide didnt work, im looking into Chop's suggestion

12/13/2008 10:05:39 PM

^^ it's possible I don't know what i'm talking about

12/13/2008 10:11:13 PM

i redfined the cond function to be the function that calculates the heat model im trying to model in 3d with my code it is:

function cond = cond(x)
c0 = 1.;
c1 = .10;
c2 = .02;
cond = c0*(1.+ c1*x + c2*x*x);

that doesnt seem to be affecting the problem cuz i was using the regular cond function just now too, yeah im still getting that
??? Error using ==> rdivide
Array dimensions must match for binary array op.

Error in ==> pcgssor at 75
alpha = rho./sum(sum(p.*q));

i tried performing an isnan divide instead but that didnt help

[Edited on December 13, 2008 at 10:29 PM. Reason : oops]

12/13/2008 10:15:57 PM

i'm pretty sure the 'error using ==> mrdivide' stuff is due to the arrays being different sizes. you need to insert some breakpoints in your pcgssor file. then you can step thru and see the different intermediate array sizes and values and how they change at each iteration.

as i see it, the problem isn't with the cond function, but the fact that the array "up" is never defined beyond an array of zeros before you perform cond(up(...)). so your ar, ag, etc arrays end up being a bunch of NaNs. when you pass them to pgssor, you end up passing a bunch of 0's and NaNs, which leads to a bunch of divide by zeros and more NaNs, and so on.

12/13/2008 10:31:39 PM

ur explaination makes sense since the method and coefficient arrays for this model are sparse
you've been a lot of help so far thanks i think ill have to wait for my professor on this last part we are all stuck on.

12/13/2008 11:20:09 PM

I think maybe the / and the ./ and the isnan./ are expecting 2d arrays and get my 3d ones, i saw somewhere with a similar function if i spilt my arrays into 2d layers i can then do the operation?

12/14/2008 5:00:47 PM