文档文库
手机版
投诉建议
热门搜索： 心得体会 演讲稿 思想汇报
首页 > 北航数值分析A大作业3

北航数值分析A大作业3

发布时间：2020-04-12   来源：文档文库   
小 中 大
字号：
手机查看
一、算法设计方案 1、解非线性方程组
将各拟合节点（xi，yj）分别带入非线性方程组，求出与(xi,yi相对应的数组te[i][j]，ue[i][j]，求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。 2、二元二次分偏插值
对数表z(t,u进行分片二次代数插值，求得对应（tij，uij）处的值，即为f(xi,yj 的值。根据给定的数表，可将整个插值区域分成 16 个小 的区域，故先判断t ij, u ij 所在，的区域，再作此区域的插值，计算 z ij，相应的Lagrange形式的插值多项式为：

p22(t,ukm1rn1l(tl°(uf(t,u
krkrm1n1其中
lk(tttwwm1tktwwkn1m1 (k=m-1, m, m+1
(r=n-1, n, n+1 l°r(uyywwn1yrywwr1


3、曲面拟合
从k=1开始逐渐增大k的值，使用最小二乘法曲面拟合法对z=f(x,y进行拟合，当107时结束计算。拟合基函数φr(xψs(y选择为φr(x=xr，ψs(y=ys。拟合系数矩阵c通过]，C(BTB1BTUG(GTG1 连续两次解线性方程组求得。C[crs其中
11B[r(xi]M1k1x0x1k1，G[s(yj]MMMMkx10Lx101x0x1LLky0y1k MMMky10Ly10y0y1LLU[f(xi,yj]
4、观察比较
计算f(xi*,y*j,p(xi*,y*j(i1,2,,8,j1,2,,5的值并输出结果，以观察p(x,y逼近f(x,y的效果。其中xi*0.1i,y*j0.50.2j。 二、全部源程序
// hean.cpp : 定义控制台应用程序的入口点。 //
#include "stdafx.h"
#include
2

#include #include
void Set_non_JacobiA(double* A,double* x;//求题中非线性方程组对应自变量向量x的雅克比//矩阵
void Set_non_B(double* B,double* x,double a,double b;//求非线性方程组Newton迭代法的右 //端式:-F(x void Array_Mult_Array(double* A,double* B,int m,int s,int n,double* C;//矩阵相乘AB＝C void Transpose(double *A,int m,int n,double* AT;//转置
void Doolittle(double *A,int n,int *M; void Solve_LU(double* A,int n,double* b,double* x;//解方程LUx＝b void Solve_lin(double* A,int n,double* B,double* x,int m;//解线性方程组Ax＝B double Vector_FanShu(double *V,int n; void Solve_non_Equation(double a,double b,double* x; //求解非线性方程组
int Near_Index(double* v,double a,int n; //查找n维向量V中与常数a最接近的元素的下标
double ChaZhi(double a,double b; //求数表z(t,u在点（x，y）处的分片二次代数差值
void NiHe(double*U,double*x,double*y,int m,int n,int p,int q,double*C; 3

//对m×n数表U(x,y进行二元多项式拟合
double Pxy(double *C,int p,int q,double x,double y; //求多项式拟合函数P(x,y在点（x，y）处的函数值

void main( { double non[4]; //非线性方程组变量向量（t,u,v,w）
double f[11][21]; //f(x,y在拟合节点处的值 double x[11],y[21]; //拟合节点 double *C; //拟合系数矩阵 double det=1e-7; //拟合精度要求 double p,d; int i,j,k,t;

//设置拟合节点
for(i=0;i<11;i++ x[i]=0.08*i; for(j=0;j<21;j++ y[j]=0.5+0.05*j;
//求拟合节点处的f(x,y的值 for(i=0;i<11;i++
{ for(j=0;j<21;j++
{ Solve_non_Equation(x[i],y[j],non; f[i][j]=ChaZhi(non[0],non[1];
} 4


} printf("\n 数值分析第三次大作业\n"; //打印数表f(xi,yi printf("\n f(xi,yi \n";

for(t=0;t<21/3;t++ {
printf("-------------------------------------------------------------------------------------\n"; printf("| f(x,y |"; for(i=3*t;i<3*t+3;i++
{ printf(" y=%4.2f |",y[i];
} printf("\n"; for(j=0;j<11;j++
{ printf("| x=%4.2f |",x[j]; for(i=3*t;i<3*t+3;i++
printf(" %+.12e |",f[j][i]; printf("\n"; }
printf("-------------------------------------------------------------------------------------\n"; printf("\n"; 5




} k=0; printf(" 不同k对应的精度 "; printf("\n----------------------------------------------"; do{ C=(double* calloc((k+1*(k+1,sizeof(double; NiHe(*f,x,y,11,21,k+1,k+1,C;

//求在当前k值拟合的节点处误差
d=0; for(i=0;i<11;i++
{ for(j=0;j<21;j++
{ p=Pxy(C,k+1,k+1,x[i],y[j]; d+=((f[i][j]-p*(f[i][j]-p; }
} printf("\n k=%d对应的σ= %.12e",k,d; if(d>=det
free(C; else {
printf("\n----------------------------------------------"; printf("\n 故可知k=k=%d<%.0e，满足题设6
要求",det; break;
} }while(++k<11;

//打印拟合系数矩阵c printf("\n\n\n 当k=%d时拟合函数p(x,y的系数矩阵c：\n",k; printf("----------------------------------------------------------------------------\n"; for(t=0;t<(k+1/3;t++
{ for(i=0;i<(k+1;i++
{ for(j=3*t;j<3*t+3;j++
printf(" %+.12e ",C[i*(k+1+j]; printf("\n"; }
printf("----------------------------------------------------------------------------\n"; }
//打印(x*,y*处的拟合逼近情况
printf("\n\n观察以下各点的逼近情况\n"; printf(" | (x*,y* | f(x*,y* | p(x*,y* | |f-p| |\n"; printf("-----------------------------------------------------------------------\n"; 7

for(i=1;i<=8;i++ for(j=1;j<=5;j++
{ Solve_non_Equation(0.1*i,0.5+0.2*j,non;
d=ChaZhi(non[0],non[1];
p=Pxy(C,k+1,k+1,0.1*i,0.5+0.2*j; printf(" | (%3.1f,%3.1f| %+.12e | %+.12e | %f |\n",0.1*i,0.5+0.2*j,d,p,abs(d-p; }
printf("-----------------------------------------------------------------------\n\n\n"; }
//求题中非线性方程组对应自变量向量x的雅克比矩阵 void Set_non_JacobiA(double* A,double* x//A为*4阶系数矩阵，x为自变量(t,u,v,w { A[0]=-0.5*sin(x[0]; A[1]=1; A[2]=1;
A[3]=1; A[4]=1; A[5]=0.5*cos(x[1]; A[6]=1;
A[7]=1; A[8]=0.5; A[9]=1; A[10]=-sin(x[2];
A[11]=1; A[12]=1; A[13]=0.5; A[14]=1; A[15]=cos(x[3]; }

8

//求非线性方程组Newton迭代法的右端式:-F(x void Set_non_B(double* B,double* x,double a,double b//a非线性方程的参数，对应题中的x //b非线性方程的参数，对应题中的y { B[0]=-(0.5*cos(x[0]+x[1]+x[2]+x[3]-a-2.67; B[1]=-(x[0]+0.5*sin(x[1]+x[2]+x[3]-b-1.07; B[2]=-(0.5*x[0]+x[1]+cos(x[2]+x[3]-a-3.74; B[3]=-(x[0]+0.5*x[1]+x[2]+sin(x[3]-b-0.79; }

//矩阵相乘AB＝C，其中A为m×s阶，B为s×n阶。
void Array_Mult_Array(double* A,double* B,int m,int s,int n,double* C { int i,j,k; for(i=0;i for(j=0;j
{
C[i*n+j]=0; for(k=0;k
C[i*n+j]+=A[i*s+k]*B[k*n+j]; } } //将m×n阶矩阵A的转置为AT void Transpose(double *A,int m,int n,double* AT { 9

int i,j; for(i=0;i for(j=0;j
AT[j*m+i]=A[i*n+j]; }

//对n阶矩阵A进行选主元的Doolittle分解
void Doolittle(double *A,int n,int *M//将分解得到的L、U仍然存储在原来的矩阵A中
{ int i,j,k,t; double *s; double Maxs,temp; s=(double* calloc(n,sizeof(double;
for(k=0;k { for(i=k;i
{
s[i]=A[i*n+k]; for(t=0;ts[i]-= A[i*n+t] * A[t*n+k]; }

Maxs=abs(s[k]; M[k]=k; for(i=k+1;i
{ if(Maxs {
10

Maxs=abs(s[i]; M[k]=i; } } if(M[k]!=k { for(t=0;t{
temp=A[k*n+t]; A[k*n+t]=A[M[k]*n+t]; A[M[k]*n+t]=temp;
} temp=s[k];
s[k]=s[M[k]]; s[M[k]]=temp;
} if(Maxs<(1e-14
{
s[k]=1e-14; printf("%.16e方阵奇异\n",Maxs; } A[k*n+k]=s[k]; for(j=k+1;(j{ for(t=0;t
A[k*n+j]-=A[k*n+t]*A[t*n+j]; A[j*n+k]=s[j]/A[k*n+k];






11
} } }

//解方程LUx＝b void Solve_LU(double* A,int n,double* b,double* x { int i,t; for(i=0;i{
x[i]=b[i]; for(t=0;tx[i]-=A[i*n+t]*x[t];
} for(i=n-1;i>-1;i--
{ for(t=i+1;t x[i]-=A[i*n+t]*x[t]; x[i]/=A[i*n+i]; } } //解线性方程组Ax＝B void Solve_lin(double* A,int n,double* B,double* x,int m//B为n×m矩阵
{ int* M,i,j; M=(int* calloc(n,sizeof(int; double *BT,*xT,temp; BT=(double* calloc(n*m,sizeof(double; 12

xT=(double* calloc(n*m,sizeof(double;

Transpose(B,n,m,BT;

Doolittle(A,n,M;//将A三角分解

for(i=0;i{ for(j=0;j
{
temp=BT[i*n+j]; BT[i*n+j]=BT[i*n+M[j]]; BT[i*n+M[j]]=temp;

}//将B转置，使得同一方程组对应的系数连续存储 Solve_LU(A,n,BT+i*n,xT+i*n;
} Transpose(xT,m,n,x; }

//求n维向量V的无穷范数
double Vector_FanShu(double *V,int n { int i; double max=0; for(i=0;i if(maxmax=abs(V[i]; return max; } 13


//求解非线性方程组
void Solve_non_Equation(double a,double b,double* x { double A[4][4],F[4],detx[4]; double det=1e-12; //求解精度要求 int i,k=0; int M=5000; //最大迭代次数 //设定迭代初始值
x[0]=0.5; x[1]=1; x[2]=1; x[3]=1; do
{ Set_non_B(F,x,a,b; Set_non_JacobiA(*A,x; Solve_lin(*A,4,F,detx,1;
if(Vector_FanShu(detx,4/Vector_FanShu(x,4
return; for(i=0;i<4;i++ x[i]+=detx[i];
k++; }while(k printf("Newton法在该初值不收敛\n"; }

//查找n维向量V中与常数a最接近的元素的下标
14

int Near_Index(double* v,double a,int n { int i,Index; double min; min=abs(v[0]-a;
Index=0; for(i=1;i{ if(min>abs(v[i]-a
{ min=abs(v[i]-a; Index=i; }
} return Index; }
//求数表z(t,u在点（x，y）处的分片二次代数差值 double ChaZhi(double a,double b { double z[6][6]={-0.5, -0.34, 0.14, 0.94, 2.06, 3.5,

-0.42, -0.5, -0.26, 0.3, 1.18, 2.38, -0.18, -0.5, -0.5, -0.18, 0.46, 1.42, 0.22, -0.34, -0.58, -0.5, -0.1, 0.62, 0.78, -0.02, -0.5, -0.66, -0.5, -0.02, 1.5, 0.46, -0.26, -0.66, -0.74,
-0.5}; double x[6]={ 0, 0.2, 0.4, 0.6, 0.8, 1 }; 15

double y[6]={ 0, 0.4, 0.8, 1.2, 1.6, 2 }; double L[3],L_[3],Z; int i,j,k,r,t; i=Near_Index(x,a,6; j=Near_Index(y,b,6; //插值区域边界处插值节点的选取 if(i==0 i=1; if(i==5 i=4; if(j==0 j=1; if(j==5 j=4; for(k=i-1;k<=i+1;k++ {
L[k-i+1]=1; for(t=i-1;t<=i+1;t++
{ if(t!=k L[k-i+1]*=((a-x[t]/(x[k]-x[t]; } }

for(r=j-1;r<=j+1;r++ {
L_[r-j+1]=1; for(t=j-1;t<=j+1;t++
{ if(t!=r L_[r-j+1]*=((b-y[t]/(y[r]-y[t];
}


16
}

Z=0; for(k=i-1;k<=i+1;k++
{ for(r=j-1;r<=j+1;r++
{ Z+=(L[k-i+1]*L_[r-j+1]*z[k][r]; }
} return Z; } //对m×n数表U(x,y进行二元多项式拟合
void NiHe(double*U,double*x,double*y,int m,int n,int p,int q,double*C //U为拟合数据点处的函数值
{ int i,j; double
*B,*BT,*BTB,*G,*GT,*GTG,*BTUG,*temp,*tempT,*CT; B=(double* calloc(m*p,sizeof(double; BT=(double* calloc(p*m,sizeof(double; BTB=(double* calloc(p*p,sizeof(double; G=(double* calloc(n*q,sizeof(double; GT=(double* calloc(q*n,sizeof(double; GTG=(double* calloc(q*q,sizeof(double; BTUG=(double* calloc(p*q,sizeof(double; temp=(double* calloc(p*n,sizeof(double; 17

for(i=0;i{ for(j=0;jB[i*p+j]=pow(x[i],j; } for(i=0;i{ for(j=0;jG[i*q+j]=pow(y[i],j; } Transpose(B,m,p,BT; Transpose(G,n,q,GT; Array_Mult_Array(BT,B,p,m,p,BTB; Array_Mult_Array(GT,G,q,n,q,GTG; Array_Mult_Array(BT,U,p,m,n,temp; Array_Mult_Array(temp,G,p,n,q,BTUG;
free(B; free(BT; free(G; free(GT; free(temp; temp=(double* calloc(p*q,sizeof(double;
Solve_lin(BTB,p,BTUG,temp,q;
free(BTB;


18
free(BTUG;
tempT=(double* calloc(q*p,sizeof(double; CT=(double* calloc(q*p,sizeof(double;
Transpose(temp,p,q,tempT; Solve_lin(GTG,q,tempT,CT,p; Transpose(CT,q,p,C;
free(temp; free(CT; free(tempT; free(GTG; }
//求多项式拟合函数P(x,y在点（x，y）处的函数值 double Pxy(double *C,int p,int q,double x,double y { int i,j; double sum=0; for(i=0;i{ for(j=0;j
sum+=(C[i*q+j]*pow(x,i*pow(y,j;
} return sum; } 19

三、输出结果




20





21



22

本文来源：https://www.2haoxitong.net/k/doc/c1bcc52a094e767f5acfa1c7aa00b52acec79c7f.html

《北航数值分析A大作业3.doc》
将本文的Word文档下载到电脑，方便收藏和打印
推荐度：
点击下载文档

文档为doc格式
相
关
案
例

相关推荐

推荐内容

【百分题库合集】消防安全知识题库 (完整word版)熵变的计算 第六章 人本主义人格流派 第1章 人体生理学绪论第三节 人体生理机能的调节 笼子里的小鸟初一作文 高等数学(170学时)课程 部队训练挑战书范文精选3篇(全文) (2021年整理)高考作文稿纸(标准) 国家公务员行测判断推理(类比推理)模拟试卷48(题后含答案及解析) 航天精神先进事迹 航天精神人物事迹范文
网站内容来自网络，如有侵权请联系我们，立即删除！ Copyright © 范文写作网京ICP备16605803号