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

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;

doubl …… 此处隐藏：4244字，全部文档内容请下载后查看。喜欢就下载吧 ……