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

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

数值分析上机作业3——求解非线性方程组

以及二元函数的插值拟合

1. 算法设计

对于全部的插值节点(xi,yj),i 0,1,...,10,j 0,1,...,20，带入非线性方程组中，用Newton迭代法解非线性方程组，得到(ti,uj),i 0,1,...,10,j 0,1,...,20。对(ti,uj)，在二维数表中进行插值，采用分片双二次插值法。插值过程中，先选择分片区域的中心节点，在数表中的列记为tt(0:5)，行记为uu(0:5)，中心节点记为(a,b)，生成向量t_temp(0:2)，

t_temp(0) (ti tt(a))(ti tt(a 1))/((tt(a 1) tt(a))(tt(a 1) tt(a 1)))， t_temp(1) (ti tt(a 1))(ti tt(a 1))/((tt(a) tt(a 1))(tt(a) tt(a 1)))， t_temp(2) (ti tt(a 1))(ti tt(a))/((tt(a 1) tt(a 1))(tt(a 1) tt(a)))，

同理，生成向量u_temp(0:2)，

u_temp(0) (uj uu(a))(uj uu(a 1))/((uu(a 1) uu(a))(uu(a 1) uu(a 1)))u_temp(1) (uj uu(a 1))(uj uu(a 1))/((uu(a) uu(a 1))(uu(a) uu(a 1))) u_temp(2) (uj uu(a 1))(uj uu(a))/((uu(a 1) uu(a 1))(uu(a 1) uu(a)))

记数表中以分片区域中心节点为中心的3×3的矩阵为T, 对于(ti,uj)插值结果为

(t_tempT)(T)(u_temp)。

在拟合p(x,y)

r,s 0

c

k

rs

xrys时，需要计算C (BTB) 1BTUG(GTG) 1，令

W BTB,M GTG，计算W,M时，根据对称性只需要计算对角线元素和对角线以上元

素即可，节省运算时间。于是WCM BUG E，用选主元的LU分解法求解WF E，再计算MC F，这里C,F只需要按行取元素进行运算即可，故不需要进行转置运算。

T

T

T

k从1到9依次增加，计算 的值，当 10 7时，得到达到精度的最小的k。

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

2. 打印输出结果

f(x,y):

4.46504018480e-001 3.24683262927e-001 2.10159686683e-001 3.40189556266e-003 -8.87358136384e-002

6.38015226510e-001 5.06611755146e-001 3.82176369277e-001 1.03043608316e-0012.64863491154e-001 1.54780200285e-001 5.19926834902e-002

8.40081395765e-001 6.99764165673e-001 5.66061442351e-001 3.19242138041e-001 2.06376192387e-001

1.05151509180e+000 9.02927430830e-001 7.60580266859e-001 4.95519756001e-001 3.73134042774e-001

1.27124675148e+000 1.11500201815e+000 9.64607727215e-001 6.82447678179e-001 5.51085208597e-001

1.49832105248e+000 1.33499863207e+000 1.17712512374e+000 8.78960023174e-001 7.39145108703e-001

1.73189274038e+000 1.56203457721e+000 1.39721691821e+000 1.08408753268e+000 9.36322772315e-001

1.97122178640e+000 1.79532959950e+000 1.62406711323e+000 1.29695464975e+000 1.14171810545e+000

k=1, sigma=3.22090897363e+000

k=2, sigma=4.65996003320e-003

k=3, sigma=1.72117537926e-004

k=4, sigma=3.30953430190e-006

k=5, sigma=2.54137773513e-008

C_rs:

2.02123036395e+000 -3.66842591519e+000 7.09246688452e-001 -4.15898479841e-001 6.74322022261e-002

3.19192646165e+000 -7.41209812962e-001 -2.69690653026e+000 -4.84601977266e-001 6.05908514261e-002

2.56706343343e-001 1.58096413206e+000 -4.65701259544e-001 1.00853116187e-001 -2.07687560816e-002

-2.68608304872e-001 -7.33963449843e-001 1.08429601112e+000 3.07285137544e-001 -4.68489486618e-002

2.16521800059e-001 -1.73026852965e-001 -8.41324310602e-002 -1.47683427939e-001 2.77711894906e-002

-5.54328606191e-002 1.40518220408e-001 -1.30388672239e-001

4.39171608118e-0016.24715198145e-0018.20347369475e-0011.02502405502e+0001.23780100674e+0001.45783058271e+000

8.48607444215e-0011.63095275395e+000-7.89186706497e-002-8.15635815143e-0012.55736987891e-0013.44966421960e-002

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

6.95959872154e-003

f(x*,y*):

0.194720 -0.183037 0.405979 -0.022516 0.634777 0.158801 0.878960 0.358651 1.136611 0.574980 1.406042 0.805941 1.685784 1.049881 1.974575 1.305335

p(x*,y*):

0.194730 -0.183042 0.405990 -0.022521 0.634787 0.158796 0.878970 0.358646 1.136620 0.574976 1.406051 0.805937 1.685791 1.049878 1.974581 1.305332

-3.30022941240e-003

-0.445498 -0.338221 -0.207366 -0.055253 0.115992 0.304429 0.508294 0.725992 -0.597567 -0.544438 -0.465358 -0.362680 -0.238568 -0.095016 0.066149 0.243254 -0.646460 -0.647361 -0.620271 -0.567565 -0.491434 -0.393902 -0.276834 -0.141950

-0.445500 -0.338224 -0.207369 -0.055255 0.115989 0.304426 0.508291 0.725989 -0.597559 -0.544430 -0.465350 -0.362671 -0.238560 -0.095009 0.066156 0.243261 -0.646446 -0.647348 -0.620257 -0.567551 -0.491421 -0.393890 -0.276822 -0.141939

对f(x,y)插值的数值情况：

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

f(xi,yi) p(xi,yi)的数值情况（z坐标为10-7）：

3. 源程序 1) 主函数

#include<stdio.h> #include<math.h>

#define length_x 11//x向量长度 #define length_y 21//y向量长度 #define x_length 8//x*向量长度 #define y_length 5//y*向量长度

extern double interpolation(double tt,double uu,double **table); //输入一点，输出一点插值

extern double row_multi_sum(double **B,int i,int j,int length); //矩阵列相乘

int main() {

double tt=0; double uu=0;

double table[6][6]={0}; double ff=0;

double xx[length_x]={0};

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

double yy[length_y]={0}; double sum=0.0;

double U[length_x][length_y]={0}; double B[length_x][10]={0}; double G[length_y][10]={0}; double W[10][10]={0}; double T[10][10]={0}; double M[10][10]={0}; double C[10][10]={0}; double F[10][length_y]={0}; double E[10][10]={0}; double XX[10][10]={0};

//double ff[length_x][length_y]={0}; //记录f阵 double pp[length_x][length_y]={0}; //记录p阵 int r=0; int s=0; int ii,jj,kk; int k_max=0;

FILE *fd;

fd = fopen("tmptmp.txt", "w"); ii=0; jj=0; kk=0;

table[0][0]=-0.5; table[0][1]=-0.34; table[0][2]=0.14; table[0][3]=0.94; table[0][4]=2.06; table[0][5]=3.5; table[1][0]=-0.42; table[1][1]=-0.5; table[1][2]=-0.26; table[1][3]=0.3; table[1][4]=1.18;

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

table[1][5]=2.38;

table[2][0]=-0.18; table[2][1]=-0.5; table[2][2]=-0.5; table[2][3]=-0.18; table[2][4]=0.46; table[2][5]=1.42;

table[3][0]=0.22; table[3][1]=-0.34; table[3][2]=-0.58; table[3][3]=-0.5; table[3][4]=-0.1; table[3][5]=0.62;

table[4][0]=0.78; table[4][1]=-0.02; table[4][2]=-0.5; table[4][3]=-0.66; table[4][4]=-0.5; table[4][5]=-0.02;

table[5][0]=1.5; table[5][1]=0.46; table[5][2]=-0.26; table[5][3]=-0.66; table[5][4]=-0.74; table[5][5]=-0.5;

for(ii=0;ii<length_x;ii++) xx[ii]=0.08*ii;

for(ii=0;ii<length_y;ii++) yy[ii]=0.5+0.05*ii;

for(ii=0;ii<length_x;ii++) //生成U阵 {

for(jj=0;jj<length_y;jj++) {

ite_equations(&tt,&uu,xx[ii],yy[jj]); //牛顿迭代解方程，输入一点（x,y）输出一点(t,u)

U[ii][jj]=(double)interpolation(tt,uu,(double **)table); //输入一点，输出一点插值

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

} }

fprintf(fd,"\nf(x,y):\n"); for(ii=0;ii<x_length;ii++) {

for(jj=0;jj<=y_length;jj++) {

fprintf(fd,"%.11e ",U[ii][jj]); }

fprintf(fd,"\n"); }

for(k_max=1;k_max<=6;k_max++) {

for(kk=0;kk<(k_max+1);kk++) //生成B阵 {

for(ii=0;ii<length_x;ii++) {

B[ii][kk]=pow(xx[ii],kk); } }

for(kk=0;kk<(k_max+1);kk++) { }

for(ii=0;ii<=k_max;ii++) //BT * B {

for(jj=0;jj<=k_max;jj++) {

W[ii][jj]=row_multi_sum((double **)B,ii,jj,length_x); } }

for(ii=0;ii<length_y;ii++) //生成G阵 {

G[ii][kk]=pow(yy[ii],kk); }

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

for(ii=0;ii<=k_max;ii++) //GT * G {

for(jj=0;jj<=k_max;jj++) {

M[ii][jj]=row_multi_sum((double **)G,ii,jj,length_y); }

}

for(ii=0;ii<=k_max;ii++) //F=B'*U for(jj=0;jj<length_y;jj++) {

F[ii][jj]=0;

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

F[ii][jj]=F[ii][jj]+B[kk][ii]*U[kk][jj]; }

for(ii=0;ii<=k_max;ii++) //E=F*G for(jj=0;jj<=k_max;jj++) {

E[ii][jj]=0;

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

E[ii][jj]=E[ii][jj]+F[ii][kk]*G[kk][jj]; }

Matrix_LU((double **)W,(double **)E,(double **)M,(double **)C,k_max,0);//(double **)XX,

for(ii=0;ii<length_x;ii++) for(jj=0;jj<length_y;jj++) {

sum=0;

for(r=0;r<=k_max;r++) for(s=0;s<=k_max;s++)

sum=sum+C[r][s]*pow(xx[ii],r)*pow(yy[jj],s); pp[ii][jj]=sum; }

sum=0;

for(ii=0;ii<length_x;ii++) for(jj=0;jj<length_y;jj++)

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

{

sum=sum+(pp[ii][jj]-U[ii][jj])*(pp[ii][jj]-U[ii][jj]); }

fprintf(fd,"\nk=%d, sigma=%.11e\n",k_max,sum);

if(sum<=1e-7) {

fprintf(fd,"\nC_rs:\n"); for(r=0;r<=k_max;r++) {

for(s=0;s<=k_max;s++) {

fprintf(fd,"%.11e ",C[r][s]); }

fprintf(fd,"\n"); }

break; } }

for(ii=0;ii<x_length;ii++) xx[ii]=0.1*(ii+1);

for(ii=0;ii<y_length;ii++) yy[ii]=0.5+0.2*(ii+1);

for(ii=0;ii<x_length;ii++) //生成f*阵 {

for(jj=0;jj<y_length;jj++) {

ite_equations(&tt,&uu,xx[ii],yy[jj]); //牛顿迭代解方程，输入一点（x,y）输出一点(t,u)

U[ii][jj]=(double)interpolation(tt,uu,(double **)table); //输入一点，输出一点插值 } }

fprintf(fd,"\nf(x*,y*):\n"); for(ii=0;ii<x_length;ii++)

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

}

for(jj=0;jj<y_length;jj++) {

fprintf(fd,"%-10f ",U[ii][jj]); }

fprintf(fd,"\n");

for(ii=0;ii<x_length;ii++) //生成p*阵 {

for(jj=0;jj<y_length;jj++) {

sum=0;

for(r=0;r<=k_max;r++) for(s=0;s<=k_max;s++)

sum=sum+C[r][s]*pow(xx[ii],r)*pow(yy[jj],s); pp[ii][jj]=sum; } }

fprintf(fd,"\np(x*,y*):\n"); for(ii=0;ii<x_length;ii++) {

for(jj=0;jj<y_length;jj++) {

fprintf(fd,"%-10f ",pp[ii][jj]); }

fprintf(fd,"\n"); }

fclose(fd); }

2) 矩阵2列相乘

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

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

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

double row_multi_sum(double **B,int i,int j,int length) //矩阵列相乘 {

double sum=0; double **TTT; int kk=0; TTT=B;

for(kk=0;kk<length;kk++) {

sum=sum+(TTT(kk,i))*(TTT(kk,j)); }

return sum; }

3) 9点分片双二次插值

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

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

double interpolation(double tt,double uu,double **table) //输入一点，输出一点插值 {

int ii=0; int jj=0;

int section_t=0; int section_u=0; double answer=0; double temp=0;

double temp_u[3]={0}; double temp_t[3]={0};

double my_t[6]={0}; double my_u[6]={0};

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

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

my_t[ii]=my_t[ii-1]+0.2; }

//判断输入点在哪个区域 //判断u用哪段插值

if ((uu>=0)&&(uu<0.6)) section_u=1;

else if((uu>=0.6)&&(uu<1.0)) section_u=2;

else if((uu>=1.0)&&(uu<1.4)) section_u=3;

else if((uu>=1.4)&&(uu<=2)) section_u=4; else ;

//判断t用哪段插值

if ((tt>=0)&&(tt<0.3)) section_t=1;

else if((tt>=0.3)&&(tt<0.5)) section_t=2;

else if((tt>=0.5)&&(tt<0.7)) section_t=3;

else if((tt>=0.7)&&(tt<=1.0)) section_t=4; else ;

temp_u[0]=(uu-my_u[section_u])*(uu-my_u[section_u+1])/((my_u[section_u-1]-my_u[section_u])*(my_u[section_u-1]-my_u[section_u+1]));

temp_u[1]=(uu-my_u[section_u-1])*(uu-my_u[section_u+1])/((my_u[section_u]-my_u[section_u-1])*(my_u[section_u]-my_u[section_u+1]));

temp_u[2]=(uu-my_u[section_u-1])*(uu-my_u[section_u])/((my_u[section_u+1]-my_u[section_u-1])*(my_u[section_u+1]-my_u[section_u]));

temp_t[0]=(tt-my_t[section_t])*(tt-my_t[section_t+1])/((my_t[section_t-1]-my_t[section_t])*(my_t[section_t-1]-my_t[section_t+1]));

temp_t[1]=(tt-my_t[section_t-1])*(tt-my_t[section_t+1])/((my_t[section_t]-my_t[section_t-1])*(my_t[section_t]-my_t[section_t+1]));

temp_t[2]=(tt-my_t[section_t-1])*(tt-my_t[section_t])/((my_t[section_t+1]-my_t[section_t-1])*(my_t[section_t+1]-my_t[section_t]));

answer=0;

for(ii=0;ii<3;ii++)

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

{

answer=answer+temp_t[ii]*(temp_u[0]*table((section_t-1+ii),section_u-1)+temp_u[1]*table((section_t-1+ii),section_u)+temp_u[2]*table((section_t-1+ii),section_u+1)); }

return answer; }

4) Newton迭代求解非线性方程组

#include<stdio.h> #include<math.h>

void ite_equations(double *tt,double *uu,double xx,double yy) //牛顿迭代解非线性方程组 {

double t=1; double u=1; double v=1; double w=1;

double d[4]={0}; //增量

double x=xx; double y=yy;

double F_Jacob[4][4]={0}; double F[4]={0};

double L[9][9]={0}; double U[9][9]={0};

double temp1=0; double temp2=0;

int ii=0; int jj=0; int kk=0;

for(ii=0;ii<=1e2;ii++) {

F_Jacob[0][0]=-0.5*sin(t);

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

F_Jacob[0][2]=1; F_Jacob[0][3]=1;

F_Jacob[1][0]=1;

F_Jacob[1][1]=0.5*cos(u); F_Jacob[1][2]=1; F_Jacob[1][3]=1; F_Jacob[2][0]=0.5; F_Jacob[2][1]=1;

F_Jacob[2][2]=-sin(v); F_Jacob[2][3]=1; F_Jacob[3][0]=1; F_Jacob[3][1]=0.5; F_Jacob[3][2]=1;

F_Jacob[3][3]=cos(w);

F[0]=0.5*cos(t)+u+v+w-x-2.67; F[1]=t+0.5*sin(u)+v+w-y-1.07; F[2]=0.5*t+u+cos(v)+w-x-3.74; F[3]=t+0.5*u+v+sin(w)-y-0.79;

LU((double **)L,(double **)U,(double **)F_Jacob,(double **)F); huidai((double **)L,(double **)U,(double **)F,(double **)d);

t=t+d[0]; u=u+d[1]; v=v+d[2]; w=w+d[3];

temp1=sqrt(t*t+u*u+v*v+w*w);

temp2=sqrt(d[0]*d[0]+d[1]*d[1]+d[2]*d[2]+d[3]*d[3]); }

if((temp2/temp1)<1e-12) {

*tt=t; *uu=u; break; }

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

5) 解线性方程组用到LU分解

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

//#define halfBW 3

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

#define 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 X(a,b) *((double *)X+(a)*(size)+b) //A 的宏定义，以方便操作二维数组

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

void LU(double **LL,double **UU,double **AA,double **bb) // //int a //需要选主元， {

double **L; double **U; double **A;

double A_temp[size][size]={0};

double **b;

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

int kk,ii,jj,pp,qq;

double s,s_temp=0;

double sum_answer=0;

double temp=0;

int s_count=0;

int n=size;

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

//////////////// 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

本文来源：https://www.bwwdw.com/article/bzyh.html