北航研究生数值分析上机作业 三 (报告+所有程序大全)(3)

//////////////// L=LL; U=UU; A=AA; b=bb;

///////////// kk=0; ii=0; jj=0; pp=0; qq=0;

////////////////////////////////// for(ii=0;ii<n;ii++) //对A_temp赋值 for(jj=0;jj<n;jj++)

A_temp[ii][jj]=A(ii,jj);

for(ii=0;ii<n;ii++) L(ii,ii)=1.0;

for (kk=0;kk<n;kk++) //总循环求出L U { sum_answer=0; s_count=kk;

for(pp=kk;pp<n;pp++)//doolitte 选主元 {

if (kk==0)

sum_answer=0; else

for(qq=0;qq<kk;qq++)

sum_answer=sum_answer+L(pp,qq)*U(qq,pp);

s=A_temp[pp][kk]-sum_answer; if (s>s_temp) {

s_temp=s; s_count=pp; }

else if (s<(-1.0)*s_temp) {

s_temp=(-1.0)*s; s_count=pp;

牛人报告,供大家参考,非常鄙视抄袭者。P.S. 北航数值分析吕老师讲的非常好,鼓励大家选这门课,能够学到很多东西。有问题联系作者!

} else ; }

for(pp=0;pp<n;pp++) //换位A_temp {

temp=A_temp[s_count][pp];

A_temp[s_count][pp]=A_temp[kk][pp]; A_temp[kk][pp]=temp; }

for(pp=0;pp<kk;pp++) //换位L {

temp=L(s_count,pp);

L(s_count,pp)=L(kk,pp); L(kk,pp)=temp; }

temp=b(s_count); b(s_count)=b(kk); b(kk)=temp;

for (ii=kk;ii<n;ii++) //先求一行U {

if(kk==0)

U(kk,ii)=A_temp[kk][ii]; else {

sum_answer=0;

for(qq=0;qq<kk;qq++)

sum_answer=sum_answer+L(kk,qq)*U(qq,ii);

U(kk,ii)=A_temp[kk][ii]-sum_answer; } }

for(ii=kk+1;ii<n;ii++) //求出一列L {

sum_answer=0;

for(qq=0;qq<kk;qq++)

sum_answer=sum_answer+L(ii,qq)*U(qq,kk);

L(ii,kk)=(A_temp[ii][kk]-sum_answer)/U(kk,kk);

牛人报告,供大家参考,非常鄙视抄袭者。P.S. 北航数值分析吕老师讲的非常好,鼓励大家选这门课,能够学到很多东西。有问题联系作者!

} } }

6) 回代求解方程组

#include<stdio.h> #include<math.h> #define size 4

#define L(a,b) *((double *)L+(a)*(size)+b) //L U 的宏定义，以方便操作二维数组

#define U(a,b) *((double *)U+(a)*(size)+b)

#define x(a) *((double *)x+a) //L U 的宏定义，以方便操作二维数组 #define y(a) *((double *)y+a) #define b(a) *((double *)b+a)

//#define xx(a) xx[a]

//A*uu=yy => L*U*uu=yy =>L*xx=yy && U*uu=xx

void huidai(double **LL,double **UU,double **bb,double **xx) //回代法求xx, LUx=b {

double **L; double **U; double **b; double **x;

double y[size]={0}; double sum_answer=0;

int ii=0; int jj=0;

L=LL; U=UU; b=bb; x=xx;

//迭代Ly=b y[0]=-b(0);

for(ii=1;ii<size;ii++)

牛人报告,供大家参考,非常鄙视抄袭者。P.S. 北航数值分析吕老师讲的非常好,鼓励大家选这门课,能够学到很多东西。有问题联系作者!

{ }

sum_answer=0;

for(jj=1;jj<ii;jj++) {

sum_answer=sum_answer+L(ii,jj)*y[jj]; }

y[ii]=-b(ii)-sum_answer;

//迭代Ux=y

x(size-1)=y[size-1]/U((size-1),(size-1));

for(ii=size-2;ii>=0;ii--) {

sum_answer=0;

for(jj=ii+1;jj<size;jj++) {

sum_answer=sum_answer+U(ii,jj)*x(jj); }

x(ii)=(y[ii]-sum_answer)/U(ii,ii); } }

7) 求解C=（B’*B）-1*(B’*U*G)*(G’*G)

#include<stdio.h> #include<math.h> #define size 10 //#define halfBW 3

#define L(a,b) L[a][b]//*((double *)L+(a)*(size)+b) //L U 的宏定义，以方便操作二维数组

#define U(a,b) U[a][b]//*((double *)U+(a)*(size)+b)

#define A(a,b) *((double *)A+(a)*(size)+b) //A 的宏定义，以方便操作二维数组 #define B(a,b) *((double *)B+(a)*(size)+b) //B 的宏定义，以方便操作二维数组 #define M(a,b) *((double *)M+(a)*(size)+b) //B 的宏定义，以方便操作二维数组 #define C(a,b) *((double *)C+(a)*(size)+b) //B 的宏定义，以方便操作二维数组 //#define X(a,b) *((double *)X+(a)*(size)+b) //A 的宏定义，以方便操作二维数组

#define X(a,b) X[a][b]

//#define b(ii) *((double *)b+ii)

void Matrix_LU(double **AA,double **BB,double **MM,double **CC,int