三对角方程组的追赶法求解方程组的解

#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>

int main()
{
char /*filename1[256],*/filename2[256];
// cout<<"Enter the path to input file:";
// cin>>filename1;
cout<<"Enter the path to output file:";
cin>>filename2;
// ifstream infile(filename1,ios::in|ios::nocreate);
ofstream outfile(filename2);
// if(!infile)
// {
// cout<<"Can not open input file:"<<filename1<<’
’;
// exit(1);
// }
int n,w,k;
// infile>>n;
n=50;
w=n+1;
double *a,*b,*d,*c,*p,*q,*x,*y,*Iy,*Ix;

a=new double[w];
d=new double[w];
c=new double[w];
b=new double[w];
p=new double[w];
q=new double[w];
x=new double[w];
y=new double[w];
Iy=new double[w];
Ix=new double[w];

//输入======================================================================================

/*

int i;
infile>>d[1]>>c[1];
// infile.seekg((n-3)*sizeof(double));
// infile.seekg((n-3)*sizeof(int),ios::cur);
for(k=0;k<n-3;k++)
infile>>a[0];
infile>>a[1];

for(k=2;k<n;k++)
{
// infile.seekg((k-2)*sizeof(int),ios::cur);
for(i=0;i<k-2;i++)
infile>>a[0];
infile>>a[k]>>d[k]>>c[k];
// infile.seekg((n-k-1)*sizeof(int),ios::cur);
for(i=0;i<n-k-1;i++)
infile>>a[0];
}
infile>>c[n];
// infile.seekg((n-3)*sizeof(int),ios::cur);
for(i=0;i<n-3;i++)
infile>>a[0];
infile>>a[n]>>d[n];
for(k=1;k<w;k++)
infile>>b[k];
*/
for(k=1;k<w;k++)
{
a[k]=1;
c[k]=1;
d[k]=4;
b[k]=k+1;
}
d[1]=2;d[n]=4;

//算法===========================================================================
if(d[1]<1e-8)
{
cout<<"Method failed!";
exit(1);
}
p[1]=d[1];
q[1]=c[1]/d[1];
for(k=2;k<n-1;k++)
{
p[k]=d[k]-a[k]*q[k-1];
if(p[k]<1e-8)
{
cout<<"Method failed!";
exit(2);
}
q[k]=c[k]/p[k];
}
p[n-1]=d[n-1]-a[n-1]*q[n-2];
if(p[n-1]<1e-8)
{
cout<<"Method failed!";
exit(3);
}
y[1]=b[1]/d[1]; Iy[1]=-a[1]/d[1];
for(k=2;k<n-1;k++)
{
y[k]=(b[k]-a[k]*y[k-1])/p[k];
Iy[k]=(-a[k]*Iy[k-1])/p[k];
}
y[n-1]=(b[n-1]-a[n-1]*y[n-2])/p[n-1];
Iy[k]=(-c[n-1]-a[k]*Iy[k-1])/p[k];
x[n-1]=y[n-1]; Ix[n-1]=Iy[n-1];
for(k=n-2;k>0;k--)
{
x[k]=y[k]-q[k]*x[k+1];
Ix[k]=Iy[k]-q[k]*Ix[k+1];
}
x[n]=(b[n]-c[n]*x[1]-a[n]*x[n-1])/(c[n]*Ix[1]+a[n]*Ix[n-1]+d[n]);
for(k=1;k<n;k++)
x[k]+=Ix[k]*x[n];
outfile<<"The roots are:(";
for(k=1;k<n;k++)
outfile<<x[k]<<’,’;
outfile<<x[n]<<");
";
delete a,d,c,p,q,x,y,Iy,Ix;
// infile.close();
outfile.close();
return 0;
}
output============

The roots are:(-6.4282,2.31175,0.18121,0.963411,0.965146,1.17601,1.33083,1.50067,1.66649,1.83338,1.99999,2.16667,2.33333,2.5,2.66667,2.83333,3,3.16667,3.33333,3.5,3.66667,3.83333,4,4.16667,4.33333,4.5,4.66667,4.83333,5,5.16667,5.33333,5.5,5.66667,5.83333,6,6.16667,6.33333,6.5,6.66666,6.83334,6.99997,7.16677,7.33293,7.5015,7.66108,7.85418,7.92219,8.45706,7.24957,12.5447);

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注