内容发布更新时间 : 2024/5/18 8:08:42星期一 下面是文章的全部内容请认真阅读。

《数值分析》计算实习题二

算法设计方案

1．主要计算步骤：

计算函数f(x,y)在拟合所需的节点处的函数值。 将各拟合节点（xi，yj）分别带入非线性方程组 0.5 cos t + u + v + w – x = 2.67 t + 0.5 sin u + v + w – y = 1.07 0.5t + u + cos v + w – x =3.74 t + 0.5u + v + sin w – y =0.79

解非线性方程组得解向量（tij，uij，vij，wij）。 对数表z(t,u)进行分片二次代数插值，求得对应（tij，uij）处的值，即为 f(xi，yj) 的值。

对上述拟合节点分别进行x，y最高次数为k（k=0,1,2,3…）次的多项式拟合。每次拟合后验证误差大小，直到满足要求。

2． 求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。 3． 对z(t，u)进行插值选择分片二次插值。

4． 拟合基函数φr(x)ψs(y)选择为φr(x)=xr，ψs(y)=ys。拟合系数矩阵c通过连续两次解线性方程组求得。

一． 源程序

#include \#include \#include \

void Doolittle(double *A,int n,int *M)

//功能说明：对n阶矩阵A进行选主元的Doolittle分解

//参数说明：A：欲进行分解的方阵，同时也是返回参数，分解后的结果 // 存储于A中 // n：方阵A的维数

// M；（返回参数）n维向量，记录选主元过程中行交换的次序 {

int i,j,k,t; double *s;

double Maxs,temp;

s=(double*) calloc(n,sizeof(double));

for(k=0;k

for(i=k;i

1

s[i]=A[i*n+k];

for(t=0;t

Maxs=abs(s[k]); M[k]=k; for(i=k+1;i

if(Maxs

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(\方阵奇异\\n\ }

A[k*n+k]=s[k];

for(j=k+1;(j

for(t=0;t

void Solve_LUEquation(double* A,int n,double* b,double* x) //功能说明：解方程LUx＝b，其中L、U共同存储在A中 //参数说明：A：经Doolittle分解后的方阵

2

// n：方阵A的维数

// b：方程组的右端向量

// x：（返回参数）方程组的解向量 {

int i,t;

for(i=0;i

x[i]=b[i];

for(t=0;t

for(i=n-1;i>-1;i--) {

for(t=i+1;t

void Transpose(double *A,int m,int n,double* AT) //功能说明：求m×n阶矩阵A的转置AT //参数说明：A：已知m×n阶矩阵 // m：A的行数 // n：A的列数

// AT：（返回参数）A的转置矩阵（n×m） {

int i,j;

for(i=0;i

for(j=0;j

void Solve_LEquation(double* A,int n,double* B,double* x,int m) //功能说明：解线性方程组Ax＝B，该函数可对系数矩阵相同 // 而右端向量不同的多个方程组同时求解。 //参数说明：A：方程组系数矩阵 // n：A的维数

// B：m组右端向量构成的n×m矩阵

// x：（返回参数）方程组的解，n×m矩阵 // m：不同右端向量的组数。 {

int* M,i,j;

M=(int*) calloc(n,sizeof(int)); double *BT,*xT,temp;

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

//求B的转置BT,使得对应一个方程组的右端系数可以连续存储

3