今天用C++摘了一段求解力学中常见的桁架问题的程序，结果在运行的时候总是无法获得想要的结果。不知道哪里出错了，特向各位大神请教。源程序和数据程序如下：
源程序：
#include
#include
#include
#include
#include
//以下定义函数
void INPUT(double X[], double Y[], int NCO[], double PROP[], double AL[],
int IB[], double REAC[]); //定义原始数据输入
void ASSEM(double X[], double Y[], int NCO[], double PROP[], double TK[][20],
double ELST[][5], double AL[]); //生成总刚度矩阵
void STIFF(int NEL, double X[], double Y[], int NCO[], double PROP[],
double ELST[][5], double AL[]); //生成单元刚度矩阵
void ELASS(int NEL, int NCO[], double TM[][20], double ELMAT[][5]); //单刚总刚传递
void BOUND(double TK[][20], double AL[], double REAC[], int IB[]);
void SLBSI(double A[][20], double B[], double D[], int N, int MS, int NRMX,
int NCMX);
void FORCE(int NCO[], double PROP[], double FORC[], double REAC[], double X[],
double Y[], double AL[]);
void OUTPT(double AL[], double FORC[], double REAC[]);
// 以下初始化全局变量
int NN,NE,NLN,NBN,N,MS;
double E,G;
// 定义输入输出数据
ifstream READ_IN; //声明一个输入流
ofstream WRITE_IO; //声明一个输出流
// 初始化程序参数
int NRMX=200;
int NCMX=20;
int NDF=2;
int NNE=2;
int NDFEL=NDF*NNE;
int main()
{
double X[100], Y[100], PROP[100], TK[200][20], AL[200], FORC[100], REAC[200]
, ELST[5][5], V[20];
int NCO[200], IB[60];
char file1[20], file2[20];
//
//打开所有文件
cin>>file1>>file2;
WRITE_IO.open(file1);
READ_IN.open(file2);
//输入数据
INPUT(X,Y,NCO,PROP,AL,IB,REAC);
//
//定义总体刚度矩阵
ASSEM(X,Y,NCO,PROP,TK,ELST,AL);
//
//引入边界条件
BOUND(TK,AL,REAC,IB);
//解总体方程
SLBSI(TK,AL,V,N,MS,NRMX,NCMX);
//计算外载荷数量
FORCE(NCO,PROP,FORC,REAC,X,Y,AL);
//
//输出
OUTPT(AL,FORC,REAC);
//
// 关闭所有文件
READ_IN.close();
WRITE_IO.close();
return 0;
}
// 子程序1：桁架静力分析子程序开始
//
///////////////////////////////////////////////////////////////////////
void INPUT(double X[], double Y[], int NCO[], double PROP[], double AL[],
int IB[], double REAC[]) //定义原始数据输入
{
//输入程序
//
int I,NUM,N1,IC[2],K,L,L1,L2,N2;
double W[3];
WRITE_IO.setf(ios::fixed);
WRITE_IO.setf(ios::showpoint);
WRITE_IO<<" "<<
"****************************************************"
<
//
//读取基本参数
READ_IN >> NN >> NE>> NLN >> NBN >> E;
WRITE_IO <<"\n\n 内部数据 \n\n"<<"结点数 :"
<
//
//按行读取结点坐标
for(I=1; I
{
READ_IN>>NUM>>X[NUM]>>Y[NUM];
}
for(I=1; I<=NN; I++)
{
WRITE_IO.precision(2);
WRITE_IO<
}
//按行读取单元连接
//单元属性
WRITE_IO
for(I=1;I
{
READ_IN>>NUM>>IC[0]>>IC[1]>>PROP[NUM];
WRITE_IO.precision(5);
WRITE_IO <
N1=NNE* (NUM-1);
NCO[N1+1]=IC[0];
NCO[N1+2]=IC[1];
}
//
//计算实际未知数目并清除载荷向量
N=NN*NDF;
for(I=1; I
{
REAC[I]=0.0;
AL[I]=0.0;
}
//读取结点载荷并按行存储
WRITE_IO
for(I=1; I
{
READ_IN>> NUM >>W[0] >>W[1];
WRITE_IO.precision(2);
WRITE_IO <
for(K=1; K
{
L=NDF*(NUM-1)+K;
AL[L]=W[K-1];
}
}
//读取边界结点数据，存储位置状态的indicatiors
//IN ARRAY IB, AND PRESCRIBED UNKNOWN VALUES IN ARRAY REAC
WRITE_IO
for(I=1; I
{
READ_IN >>NUM>> IC[0]>> IC[1]>> W[0]>> W[1];
WRITE_IO.precision(4);
WRITE_IO <<setw(10)<< NUM
<<setw(10)<< IC[0]
<<setw(10)<< IC[1]
<<setw(20)<< W[0]
<<setw(10)<< W[1]
<<"\n";
L1=(NDF+1)*(I-1)+1;
L2=NDF*(NUM-1);
IB[L1]=NUM;
for(K=1; K<=NDF; K++)
{
N1=L1+K;
N2=L2+K;
IB[N1]=IC[K-1];
REAC[N2]=W[K-1];
}
}
//子程序结束
return;
}
CSDN有限制，所以完整程序在http://emuch.net/bbs/viewthread.php?tid=9072399
可读数据：将以下数据新建一个*.DAT文件并保存到与程序一致的文件夹中，*可以是任意名称，之后输出用*.RES，*为同一个文件名。
6,10,2,2,2.0E8,
1,2.0,1.0
2,2.0,0.0
3,1.0,1.0
4,1.0,0.0
5,0.0,1.0
6,0.0,0.0
1,1,2,1.000
2,1,3,1.000
3,2,4,1.000
4,1,4,1.414
5,2,3,1.414
6,3,4,1.000
7,3,5,1.000
8,3,6,1.414
9,4,5,1.414
10,4,6,1.000
2,0.0,-1.0
4,0.0,-1.0
5,0,1,0.0,0.0
6,0,1,0.0,0.0
如果可以运行，那么首先输入*.RES作为结果的放置空间，再读入*.DAT文件。
不胜感激
