Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
how solve a tridiagonal equation system when we are in three dimensional space?

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: ghasem

Date: 13 Apr, 2013 10:59:07

Message: 1 of 17

hi to all.

I have one question.

I want to solve tridiagonal equation system,when a,b,c arrays in coefficient matrix of A are 3D.i,e as following:

Ax = f

or

a(i,j,k)*X(i,j-1,k)+b(i,j,k)*X(i,j,k)+c(i,j,k)*X(i,j+1,k) = f(i,j,k)

that i=1:50 , j = 1:70 , k=1:110;

also,I must mentioned that :

1- a = c

2- 3D a,b,c arrays are created and saved(as 3d full matrices with 50*70*110 size) before the equation solving,and the right hand side(f(i,j,k)) change each time that I want to solve the equation.

what method is best choice for solving?

I know that mldevide(\) method is better(faster) than Thomas algorithms.

but I want to use from sparse matrix concept for 3D case.in 2D I solved Ax=b very simple by spars matrix. but I don't know that what should I do in 3D case.

I saw "ndSparse" method for ND sparse marices,but I don't know how work with it.

I can solve this tridiagonal system with Thomas algorithm,but can not with mldevide method.

in fact,I don't know that how construct A sparse matrix from 3D a,b,c arrays (that are sub-diagonal,main-diagonal and super-diagonal respectively),and also I don't know how implement x = A\f in 3D case.

is there anyone that help me and direct me a little?

best regards...
ghasem...

Subject: how solve a tridiagonal equation system when we are in three

From: Nasser M. Abbasi

Date: 13 Apr, 2013 11:31:58

Message: 2 of 17

see comments here under this submission:

http://www.mathworks.com/matlabcentral/fileexchange/1359

> in fact,I don't know that how construct A sparse matrix from
>3D a,b,c arrays (that are sub-diagonal,main-diagonal and super-diagonal
> respectively),and also I don't knowhow implement x = A\f in 3D case.
>
> is there anyone that help me and direct me a little?
>
> best regards...
> ghasem...
>

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: Bruno Luong

Date: 13 Apr, 2013 11:56:08

Message: 3 of 17

"ghasem " <shaban_sadeghi@yahoo.com> wrote in message <kkbdpr$63p$1@newscl01ah.mathworks.com>...
> hi to all.
>
> I have one question.
>
> I want to solve tridiagonal equation system,when a,b,c arrays in coefficient matrix of A are 3D.i,e as following:
>
> Ax = f
>
> or
>
> a(i,j,k)*X(i,j-1,k)+b(i,j,k)*X(i,j,k)+c(i,j,k)*X(i,j+1,k) = f(i,j,k)
>

Your equation involves only the second dimension (j). The two other dimensions (i,k) can be considered as parameters.

In fact you have 50*110 linear systems of (70 x 70) to solve.

Of course you can also consider one big linear system of 50*110 of tridiagonal blocks (385000 x 385000). Modern computer can deal easily with such size of sparse matrix.

To build your sparse matrix, simply number your equations as row and your unknown as column. It doesn't matter where as these equations or unknowns represent a data physically in 1, 2, 3, or 100 dimensional space.

Bruno

Subject: how solve a tridiagonal equation system when we are in three

From: ghasem

Date: 13 Apr, 2013 12:09:08

Message: 4 of 17

"Nasser M. Abbasi" wrote in message <kkbfnb$7v3$1@speranza.aioe.org>...
> see comments here under this submission:
>
> http://www.mathworks.com/matlabcentral/fileexchange/1359
>
================
thank you Mr.Abbasi.
but this submission is very slow!
I want to solve my problem by fastest possible way...
I think that mldevide is good choice,but I don't konw how write x=A\b in three dimensional space.
do you any idea?

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: ghasem

Date: 13 Apr, 2013 12:29:08

Message: 5 of 17


> Your equation involves only the second dimension (j). The two other dimensions (i,k) can be considered as parameters.
>
> In fact you have 50*110 linear systems of (70 x 70) to solve.
>
> Of course you can also consider one big linear system of 50*110 of tridiagonal blocks (385000 x 385000). Modern computer can deal easily with such size of sparse matrix.
>
> To build your sparse matrix, simply number your equations as row and your unknown as column. It doesn't matter where as these equations or unknowns represent a data physically in 1, 2, 3, or 100 dimensional space.
>
> Bruno
==================
thank you very much fro your attention.
I know that I have 50*110 linear systems of (70 x 70) to solve.
but I don't now how solve it with mldevide.
for example,one section from my code is as following:
first step:
I calculated 3d full matrix of a(i,j,k) and b(i,j,k).
second step:
I want to create tridiagonal coefficient matrix (A in Ax=b).so,I wrote this code:
for i =1:nx
    for k =2:nz
        Amd = sparse(1:ny+1 ,1:ny+1 ,b(i,:,k) ,ny+1,ny+1);
        Asub = sparse(2:ny+1 ,1:ny ,a(i,2:ny+1,k),ny+1,ny+1);
        A{i,k}= Amd + Asub + Asub.' ;
    end
end
third step:
in my main program,I want to find unknowns(X).so I wrote:
X = cell2mat(A) \ f ;
==============
But in third step,I get this error:
??? Error using ==> ldivide
Number of array dimensions must match for binary array op.
in fact,I don't know how insert index in third step,truly.can you give me a little explanation?
also,I have another question.in second step,my code is good?or is there any better way?
best regards...
ghasem

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: Bruno Luong

Date: 13 Apr, 2013 12:56:07

Message: 6 of 17

"ghasem " <shaban_sadeghi@yahoo.com> wrote in message <kkbj2k$j6o$1@newscl01ah.mathworks.com>...
>
> > Your equation involves only the second dimension (j). The two other dimensions (i,k) can be considered as parameters.
> >
> > In fact you have 50*110 linear systems of (70 x 70) to solve.
> >
> > Of course you can also consider one big linear system of 50*110 of tridiagonal blocks (385000 x 385000). Modern computer can deal easily with such size of sparse matrix.
> >
> > To build your sparse matrix, simply number your equations as row and your unknown as column. It doesn't matter where as these equations or unknowns represent a data physically in 1, 2, 3, or 100 dimensional space.
> >
> > Bruno
> ==================
> thank you very much fro your attention.
> I know that I have 50*110 linear systems of (70 x 70) to solve.
> but I don't now how solve it with mldevide.
> for example,one section from my code is as following:
> first step:
> I calculated 3d full matrix of a(i,j,k) and b(i,j,k).
> second step:
> I want to create tridiagonal coefficient matrix (A in Ax=b).so,I wrote this code:
> for i =1:nx
> for k =2:nz
> Amd = sparse(1:ny+1 ,1:ny+1 ,b(i,:,k) ,ny+1,ny+1);
> Asub = sparse(2:ny+1 ,1:ny ,a(i,2:ny+1,k),ny+1,ny+1);
> A{i,k}= Amd + Asub + Asub.' ;
> end
> end
> third step:
> in my main program,I want to find unknowns(X).so I wrote:
> X = cell2mat(A) \ f ;
> ==============

Why don't you solve each linear system within the loop?

IMO it is also faster to create the sparse matrix with a single sparse command (using the indices of all tridiagonal together).

Bruno

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: ghasem

Date: 13 Apr, 2013 13:42:08

Message: 7 of 17

> Why don't you solve each linear system within the loop?
>
> IMO it is also faster to create the sparse matrix with a single sparse command (using the indices of all tridiagonal together).
>
> Bruno
===================
I think that by using "for-loop" into Ax=f systems , my naswer become very slow.
in fact,I have to solve Ax=f in all of the three directions of x,y,z.and also,these equations are within one "for-loop" for time.i.e I have:
for n=1:10000
% here is early works of my program.(I don't mentioned them here)
% next,I have to solve AU=b that U=X,Y,Z . i,e:
1-
ax(i,j,k)* X(i,j-1,k) + bx(i,j,k)* X(i,j,k) + cx(i,j,k) * X(i,j+1,k) = fx(i,j,k);
2-
ay(i,j,k)* Y(i,j,k-1) + by(i,j,k)* Y(i,j,k) + cy(i,j,k) * Y(i,j,k+1) = fy(i,j,k);
3-
az(i,j,k)* Z(i-1,j,k) + bz(i,j,k)* Z(i,j,k) + cz(i,j,k) * Z(i+1,j,k) = fz(i,j,k);
that:
i=1:50 ; j =1:70 ; k=1:110;
% here is my other sections of my program.(I don't mentioned them here)
end
how solve these equations,so that the solution to be best method?
===================
also,I must to mentioned that when I used from thomas algorithm method,my program took 80 minutes!while when I solved my program by other method(i,e without Ax=b systems) it took only 8 minutes!
and I have to use from Ax=b method in my program,also(because this is a research).
best regards...
ghasem...

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: Bruno Luong

Date: 13 Apr, 2013 14:02:08

Message: 8 of 17

"ghasem " <shaban_sadeghi@yahoo.com> wrote in message <kkbnbg$he$1@newscl01ah.mathworks.com>...
> > Why don't you solve each linear system within the loop?
> >
> > IMO it is also faster to create the sparse matrix with a single sparse command (using the indices of all tridiagonal together).
> >
> > Bruno
> ===================
> I think that by using "for-loop" into Ax=f systems , my naswer become very slow.

Nonsense, it take 0.24 second on my old laptop.

ny = 70;
nxnz = 50*110;
a = rand(ny-1,nxnz);
b = rand(ny,nxnz);
c = rand(ny-1,nxnz);
i = [2:ny 1:ny 1:ny-1]';
j = [1:ny-1 1:ny 2:ny]';
s = [a; b; c];
f = rand(ny,nxnz);
X = zeros(ny,nxnz);
tic
for k=1:nxnz
    S = sparse(i(:), j(:), s(:,k), ny, ny);
    X(:,k) = S \ f(:,k);
end
t = toc

% Bruno

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: ghasem

Date: 13 Apr, 2013 14:25:08

Message: 9 of 17

> Nonsense, it take 0.24 second on my old laptop.
>
> ny = 70;
> nxnz = 50*110;
> a = rand(ny-1,nxnz);
> b = rand(ny,nxnz);
> c = rand(ny-1,nxnz);
> i = [2:ny 1:ny 1:ny-1]';
> j = [1:ny-1 1:ny 2:ny]';
> s = [a; b; c];
> f = rand(ny,nxnz);
> X = zeros(ny,nxnz);
> tic
> for k=1:nxnz
> S = sparse(i(:), j(:), s(:,k), ny, ny);
> X(:,k) = S \ f(:,k);
> end
> t = toc
>
> % Bruno
============================
thank you very much.you code is very nice!
I'll try with your pattern,now and I'll declare my answer,as soon as.
ghasem

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: ghasem

Date: 13 Apr, 2013 15:01:16

Message: 10 of 17

> Nonsense, it take 0.24 second on my old laptop.
>
> ny = 70;
> nxnz = 50*110;
> a = rand(ny-1,nxnz);
> b = rand(ny,nxnz);
> c = rand(ny-1,nxnz);
> i = [2:ny 1:ny 1:ny-1]';
> j = [1:ny-1 1:ny 2:ny]';
> s = [a; b; c];
> f = rand(ny,nxnz);
> X = zeros(ny,nxnz);
> tic
> for k=1:nxnz
> S = sparse(i(:), j(:), s(:,k), ny, ny);
> X(:,k) = S \ f(:,k);
> end
> t = toc
>
> % Bruno
====================
sorry Bruno.I have another question.
in your code,How convert X from sparse to full matrix,again?
Since,in next time loop of program,I have to re-calculate RHS of equation and some other unknowns by new values of X.
the only method for my purpose, is "full" command?
thank you very much for your attention...
ghasem...

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: Bruno Luong

Date: 13 Apr, 2013 16:02:08

Message: 11 of 17

"ghasem " <shaban_sadeghi@yahoo.com> wrote in message


> in your code,How convert X from sparse to full matrix,again?
>

> > X = zeros(ny,nxnz);

In my code X preallocated with ZEROS() is a _full_ matrix. It never becomes sparse.

Bruno

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: ghasem

Date: 13 Apr, 2013 17:10:08

Message: 12 of 17

> > in your code,How convert X from sparse to full matrix,again?
> >
>
> > > X = zeros(ny,nxnz);
>
> In my code X preallocated with ZEROS() is a _full_ matrix. It never becomes sparse.
>
> Bruno
=======================
I know that X preallocate with ZEROS() .But I think that,when we calculate
 X(:,k) = S \ f(:,k)
in fact first dimension of X is j and original X(i,j,k) is replaced with X(j,i*k).
do I think truly?
on the other hand,
other calculations in program that need to X(i,j,k),don't disrupted?
for example,after tridiagonal system to be solved,I have to calculate:
H(i,j,k) = q(i,j,k)*(X(i,j,k) - X(i,j-1,k))+...
that q(i,j,k) coefficient is determined.
So,H(i,j,k) values are correct and I don't need to reshape X(:,1:nxnz) values obtained from tridiagonal system solving,to X(i,j,k)?
my question is clear?
best wishes...
ghasem

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: ghasem

Date: 14 Apr, 2013 06:27:08

Message: 13 of 17

> I know that X preallocate with ZEROS() .But I think that,when we calculate
> X(:,k) = S \ f(:,k)
> in fact first dimension of X is j and original X(i,j,k) is replaced with X(j,i*k).
> do I think truly?
> on the other hand,
> other calculations in program that need to X(i,j,k),don't disrupted?
> for example,after tridiagonal system to be solved,I have to calculate:
> H(i,j,k) = q(i,j,k)*(X(i,j,k) - X(i,j-1,k))+...
> that q(i,j,k) coefficient is determined.
> So,H(i,j,k) values are correct and I don't need to reshape X(:,1:nxnz) values obtained from tridiagonal system solving,to X(i,j,k)?
> my question is clear?
> best wishes...
> ghasem
========================================
please help me...
when in my code,I used from mldevide and sparse matrix to find unknown values,(say X,Y,Z),the other unknowns in my program that are depend on X,Y,Z values,are not changed.for example,I have to calculate:
H = coeff1* (Y(1:nx+1,1:ny,2:nz+1) - Y(1:nx+1,1:ny,1:nz)) ...
    + coeff2* (Z(1:nx+1,2:ny+1,1:nz) - Z(1:nx+1,1:ny,1:nz));

But when I used from Bruno method(in above),H value didn't update.
I think that when I calculate X(:,k)( or Y or Z),in fact size of original X(i,j,k),Y(i,j,k),Z(i,j,k) is changed and they become 2D (no 3D).
How re-convert X(:,k) to X(i,j,k)to continue ordinary operations in continue?
best regards...
ghasem

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: Bruno Luong

Date: 14 Apr, 2013 08:06:11

Message: 14 of 17

"ghasem " <shaban_sadeghi@yahoo.com> wrote in message <kkdi7s$166$1@newscl01ah.mathworks.com>...
> > I know that X preallocate with ZEROS() .But I think that,when we calculate
> > X(:,k) = S \ f(:,k)
> > in fact first dimension of X is j and original X(i,j,k) is replaced with X(j,i*k).
> > do I think truly?
> > on the other hand,
> > other calculations in program that need to X(i,j,k),don't disrupted?
> > for example,after tridiagonal system to be solved,I have to calculate:
> > H(i,j,k) = q(i,j,k)*(X(i,j,k) - X(i,j-1,k))+...
> > that q(i,j,k) coefficient is determined.
> > So,H(i,j,k) values are correct and I don't need to reshape X(:,1:nxnz) values obtained from tridiagonal system solving,to X(i,j,k)?
> > my question is clear?
> > best wishes...
> > ghasem
> ========================================
> please help me...
> when in my code,I used from mldevide and sparse matrix to find unknown values,(say X,Y,Z),the other unknowns in my program that are depend on X,Y,Z values,are not changed.for example,I have to calculate:
> H = coeff1* (Y(1:nx+1,1:ny,2:nz+1) - Y(1:nx+1,1:ny,1:nz)) ...
> + coeff2* (Z(1:nx+1,2:ny+1,1:nz) - Z(1:nx+1,1:ny,1:nz));
>
> But when I used from Bruno method(in above),H value didn't update.
> I think that when I calculate X(:,k)( or Y or Z),in fact size of original X(i,j,k),Y(i,j,k),Z(i,j,k) is changed and they become 2D (no 3D).
> How re-convert X(:,k) to X(i,j,k)to continue ordinary operations in continue?

help reshape
help permute

Bruno

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: ghasem

Date: 14 Apr, 2013 09:39:08

Message: 15 of 17

> help reshape
> help permute
>
> Bruno
=====================
I know that I can use from reshape and permute.but because my program is within a time-loop (for n=1:10000),so when I use from reshape and permute,repeatedly,my program is not fast enough.and I must use several variable and preallocate them.So,this method take a lot of memory.
there is not any other method that I can solve tridiagonal equation system,without use from permute and reshape,repeatedly and also to be fast?
i,e mldevide and reshape and permute are best and fastest method?
best regards
ghasem

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: Bruno Luong

Date: 14 Apr, 2013 09:55:07

Message: 16 of 17

"ghasem " <shaban_sadeghi@yahoo.com> wrote in message <kkdtfs$stp$1@newscl01ah.mathworks.com>...
> > help reshape
> > help permute
> >
> > Bruno
> =====================
> I know that I can use from reshape and permute.but because my program is within a time-loop (for n=1:10000),so when I use from reshape and permute,repeatedly,my program is not fast enough.

Why you should do it repeatly?

>> A = rand(50,70,200);
>>tic; B = permute(A,[2 1 3]); B = reshape(B, 70, []); toc
Elapsed time is 0.007224 seconds.

It's 60000 faster than 8 minutes. A peanut in the ocean.

Please post the small example running code when asking the question. I can't afford to constantly guess what exactly _you_ do.

Bruno

Subject: how solve a tridiagonal equation system when we are in three dimensional space?

From: ghasem

Date: 14 Apr, 2013 11:35:08

Message: 17 of 17

> Why you should do it repeatly?
>
> >> A = rand(50,70,200);
> >>tic; B = permute(A,[2 1 3]); B = reshape(B, 70, []); toc
> Elapsed time is 0.007224 seconds.
>
> It's 60000 faster than 8 minutes. A peanut in the ocean.
>
> Please post the small example running code when asking the question. I can't afford to constantly guess what exactly _you_ do.
>
> Bruno
==================================
sorry Bruno,you are right.I write small example fro my code:
first step:
I calculate: aq(i,j,k),bq(i,j,k),cq(i,j,k), that q=u,v,w,x,y,z and i=1:50,j=1:70,k=1:110.
second step:
preallocate U,V,W,X,Y,Z =>ie : = zeros(50,70,110)
third step:
I have a time loop that within this loop,I have to solve Ax=b for U,V,W,X,Y,Z unknowns.
and also,RHS of this six systems of Ax=b,are depend on new values of two other unknowns.
==============================
i.e,one section from my code is as following:
for time =1:10000

% for X:
fx(1:nx,2:ny,2:nz) = coeff1* (W(1:nx,2:ny,2:nz) - W (1:nx,1:ny-1,2:nz))...
                            +coeff2* (V(1:nx,2:ny,2:nz) - V (1:nx,2:ny,1:nz-1));

fx2 = permute(fx ,[2,1,3]);
fx3 = reshape(fx2,[ny+1,nx*(nz+1)]);
for k=1:nxnze
    Ax2 = sparse(i(:), j(:), Ax(:,k), ny+1, ny+1);
    X1(:,k) = Ax2 \ fX3(:,k);
end
X3 = reshape(X2,[ny+1,nx,(nz+1)]);
X = permute(X3,[2 1 3]);
==============================
this was Ax=f for X unknown.
now,suppose I have five other tridiagonal system.

forth step:
% now after Ax=f for x =U,V,W,X,Y,Z , I have to calculate smaller section,aslo.
% for example,for X,I have to calculate as following,again:

fx4(:,1:10,:) = coeff3*(U(:,2:11,:) - U(:,1:10,:));
fx5 = permute(fx4 ,[2,1,3]);
fx6 = reshape(fx5 ,[10 ,nx*(nz+1)]);
for k=1:nxnze
    Ax4 = sparse(ii(:), jj(:) , Ax3(:,k) ,10 ,10);
    X4 (:,k) = A_ex4 \ fx6(:,k);
end
X5 = reshape(X4 ,[10 , nx , (nz+1)]);
X6 = permute(X5 ,[2 1 3]);
%
X(:, 2:11 ,:) = X(: , 2:11 ,:) + X6 ;

end % related to for time=...

% NOW suppose I have same procedure for U,V,W,Y,Z.and all of these calculations are within for-loop of time.So,I want to use from best method for Ax=f.such that my program to be fast.

what is your idea that my program to be fast?

Mr,Bruno,can I have your email?I want to tell you that what is my purpose.
best regard...
ghasem...

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us