solving large system of equations

vsgdp · science forum addict Joined: 01 May 2005 Posts: 78

Hi,

This is a follow up to my "help with finite difference scheme" post.
Instead of doing my hack to create an explicit scheme, I just left it in
implicit form and solve a system of equations at every step (which is
working now!). Right now I am taking the inverse to solve the system.
However, it is taking some time. I know there are better methods for
solving. My system is tridiagonal except for the first and last row of the
matrix. The system looks like this.

[1 c 0 ... 0 -c ]
..
..
..
[0 ... 0 -c 1 c 0 ... 0]
..
..
..
[c 0 0 ... -c 1]

In particular, I am looking for a MATLAB function that will solve this
efficiently. I don't really want to code up my own right now since the
course topic is on numeric differential equations and not matrix
computations.

Ron Shepard · science forum beginner Joined: 11 Jun 2005 Posts: 31

In article <wDEqe.6804$6s.4099@fed1read02>, "vsgdp" <spam@null.com>
wrote:

spasmous · science forum addict Joined: 03 May 2005 Posts: 66

MATLAB has lsqr() implemented, if you don't need anything particularly
specialized - ie. not optimal for your matrix.

Carl Barron · science forum beginner Joined: 12 Jun 2005 Posts: 27

Ron Shepard <ron-shepard@NOSPAM.comcast.net> wrote:

spasmous · science forum addict Joined: 03 May 2005 Posts: 66

Carl Barron wrote:

Peter Spellucci · science forum Guru Joined: 29 Apr 2005 Posts: 702

In article <wDEqe.6804$6s.4099@fed1read02>,
"vsgdp" <spam@null.com> writes:

Jean-Claude Arbaut · science forum Guru Joined: 13 Jun 2005 Posts: 573

Some C code I needed for the same job (finite difference). Comes with
no guarantee :-)

/*
tridiag(n,a,b,c,d)

b[i] * x[i-1] + a[i] * x[i] + c[i] * x[i+1] = d[i], 0 <= i < n
b[0] = c[n-1] = 0

Solution in a[0] ... a[n-1]
*/

void tridiag(int n,double *a,double *b,double *c,double *d) {
int i;
double t;

for(i=1;i<n;i++) {
t=b[i]/a[i-1];
a[i]-=c[i-1]*t;
d[i]-=d[i-1]*t;
}
a[n-1]=d[n-1]/a[n-1];
for(i=n-2;i>=0;i--) a[i]=(d[i]-c[i]*a[i+1])/a[i];
}

/*
tricyc(n,a,b,c,d,u,v)

b[i] * x[i-1] + a[i] * x[i] + c[i] * x[i+1] + u[i] * x[n-1] = d[i], 0 <=
i < n

add to the first row: b[0] * x[n-1]
add to the last row: c[n-1] * x[0] + v[0] * x[0] + ... + v[n-1] * x[n-1]

Solution in a[0] ... a[n-1]
*/

void tricyc(int n,double *a,double *b,double *c,double *d,double *u,double
*v) {
int i;
double t,z;

u[0]+=b[0];
v[0]+=c[n-1];
for(i=1;i<n;i++) {
t=b[i]/a[i-1];
a[i]-=c[i-1]*t;
u[i]-=u[i-1]*t;
d[i]-=d[i-1]*t;
z=v[i-1]/a[i-1];
v[i]-=c[i-1]*z;
v[n-1]-=u[i-1]*z;
d[n-1]-=d[i-1]*z;
}
a[n-1]=d[n-1]/(a[n-1]+u[n-1]+v[n-1]);
for(i=n-2;i>=0;i--) a[i]=(d[i]-c[i]*a[i+1]-u[i]*a[n-1])/a[i];
}

/*
penta(n,a,b,c,d,e,z)

a c e .
b \ \ \
d \ \ \
. \ \ \

Solution in a
*/

void penta(int n,double *a,double *b,double *c,double *d,double *e,double
*z) {
int i;
double t,u;

for(i=2;i<n;i++) {
t=b[i-1]/a[i-2];
a[i-1]-=c[i-2]*t;
c[i-1]-=e[i-2]*t;
z[i-1]-=z[i-2]*t;
u=d[i]/a[i-2];
b[i]-=c[i-2]*u;
a[i]-=e[i-2]*u;
z[i]-=z[i-2]*u;
}
t=b[n-1]/a[n-2];
a[n-1]-=c[n-2]*t;
z[n-1]-=z[n-2]*t;
a[n-1]=z[n-1]/a[n-1];
a[n-2]=(z[n-2]-c[n-2]*a[n-1])/a[n-2];
for(i=n-3;i>=0;i--) a[i]=(z[i]-c[i]*a[i+1]-e[i]*a[i+2])/a[i];
}

Le 13/06/2005 15:38, dans d8k285$lng$1@fb04373.mathematik.tu-darmstadt.de,
« Peter Spellucci » <spellucci@fb04373.mathematik.tu-darmstadt.de> a écrit :

Jean-Claude Arbaut · science forum Guru Joined: 13 Jun 2005 Posts: 573

Ok, the answer was still truncated... Sorry Smile

I hope this time it's OK.

/*
penta(n,a,b,c,d,e,z)

a c e *
b \ \ \
d \ \ \
* \ \ \

Solution in a
*/

void penta(int n,double *a,double *b,double *c,double *d,double *e,double
*z) {
int i;
double t,u;

for(i=2;i<n;i++) {
t=b[i-1]/a[i-2];
a[i-1]-=c[i-2]*t;
c[i-1]-=e[i-2]*t;
z[i-1]-=z[i-2]*t;
u=d[i]/a[i-2];
b[i]-=c[i-2]*u;
a[i]-=e[i-2]*u;
z[i]-=z[i-2]*u;
}
t=b[n-1]/a[n-2];
a[n-1]-=c[n-2]*t;
z[n-1]-=z[n-2]*t;
a[n-1]=z[n-1]/a[n-1];
a[n-2]=(z[n-2]-c[n-2]*a[n-1])/a[n-2];
for(i=n-3;i>=0;i--) a[i]=(z[i]-c[i]*a[i+1]-e[i]*a[i+2])/a[i];
}

Google