square.cpp

数值计算方法中几个重要的算法用VC++实现

/*--解对称正定矩阵方程组的平方根法--*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
double *Assign(int n)
{
	return (double*)malloc(sizeof(double)*(n+1));
}
void main()
{
	int i,ik,jk,j,k,n;/*--n是未知数的个数--*/
    double **A,*b,*x,*y,*l,sum;
    printf("Input dimension of Array n(n>=1):\n");
	scanf("%d",&n);
	/*--为数组A,b,x,y,l--*/
    A=(double **)malloc(sizeof(double*)*(n+1));
	for(i=1;i<=n;i++)
		A[i]=Assign(n);
	b=Assign(n);x=Assign(n);y=Assign(n);l=Assign(n);
	/*--对数组A,b进行赋值--*/
	printf("Input the element of Array A:\n");
	for(i=1;i<=n;i++)
		for(j=1;j<=n;j++)
			scanf("%lf",&A[i][j]);
	printf("Input the element of Array b:\n");
	for(i=1;i<=n;i++)
		scanf("%lf",&b[i]);
	/*--A=LL(t)--分解计算--*/
	l[1]=sqrt(A[1][1]);
	for(i=2;i<=n;i++)
	{
		jk=i*(i-1)/2+1;
		l[jk]=A[i][1]/l[1];
	}
	/*--对于j=2,3,……,n--*/
	for(j=2;j<=n;j++)
	  {
	  	sum=0;
		for(k=1;k<=j-1;k++)
		{
			jk=j*(j-1)/2+k;
			sum+=pow(l[jk],2);
		}
		l[j*(j-1)/2+j]=sqrt(A[j][j]-sum);
		if(j!=n)
		{
			for(i=j+1;i<=n;i++)
			{
				sum=0;
				for(k=1;k<=j-1;k++)
				{
					ik=i*(i-1)/2+k;
					jk=j*(j-1)/2+k;
					sum+=l[ik]*l[jk];
				}
				l[i*(i-1)/2+j]=(A[i][j]-sum)/l[j*(j-1)/2+j];
			}
		}/*-if-*/
	}/*-for(j)-*/
	/*--求解Ly=b--*/
	y[1]=b[1]/l[1];
	for(i=2;i<=n;i++)
	{
		sum=0;
		for(k=1;k<=i-1;k++)
		{
			ik=i*(i-1)/2+k;
			sum+=l[ik]*y[k];
		}
		y[i]=(b[i]-sum)/l[i*(i-1)/2+i];
	}
	/*--求解L(t)x=y--*/
	x[n]=y[n]/l[n*(n-1)/2+n];
	for(i=n-1;i>=1;i--)
	{
		sum=0;
		for(k=i+1;k<=n;k++)
		{
			jk=k*(k-1)/2+i;
			sum+=l[jk]*x[k];
		}
		x[i]=(y[i]-sum)/l[i*(i-1)/2+i];
	}
	/*-输出x1,x2,x3.....,xn-*/
	for(i=1;i<=n;i++)
		printf("x%d=%.8lf\n",i,x[i]);
	/*--释放空间--*/
    free(l);free(y);free(x);free(b);
	for(i=n;i>=1;i--)
		free(A[i]);
	
}