北航数值分析第三次大作业

数值分析第三次大作业

一、算法的设计方案 1、求解非线性方程组

将题目中给出的(,)i i x y 当作已知量代入题目给定的非线性方程组，求出与

(,)i i x y 相对应的数组te[i][j]，ue[i][j]，此处采用的是牛顿法解非线性方程组，其

算法如书上91页所示。 2、分片二次代数插值

对所求出的数组te[i][j]，ue[i][j]，通过分片二次代数插值运算，得到与数组te[11][21]，ue[11][21]对应的数组ze[11][21]，从而得到二元函数z=(,)i i f x y ，此处采用如书上101页例2中所示的分片二次代数插值。 3、曲面插值

利用x[11]，y[21]，ze[11][21]建立二维函数表，进行曲面插值计算，逐步提高k 值，计算其精度，看其是否满足要求，以此来确定循环结束的时刻，并得到曲面拟合的系数矩阵C[r][s]，此处的算法如书142页所示，只需将所需矩阵给出，然后按公式进行计算即可。 4、比较

观察和),(j i y x p 逼近(,)i i f x y 的效果。观察逼近效果只需要利用新给的点列

(,)i i x y 重复上面（1）和（2）的过程，得到与新的插值节点(,)i i x y 对应的(,)i i f x y ，

再与对应的(,)i i p x y 比较即可，这里求解(,)i i p x y 可以直接使用（3）中的C[r][s]和k 。 5、几点说明

分片二次插值的结果x[i]，y[j]，ze[i][j]输出到一个文件shubiao.txt 中，方便结果的复制与粘贴。

曲面插值的结果输出到一个文件xishu.txt 中，包括循环中每一次的k 值以及误差平方和sigma 的值，还有最后满足误差要求时曲面插值的系数C[r][s]。

观察逼近效果的结果输出到一个文件shubiao1.txt 中，方便结果的复制与粘贴。

二、源程序如下：

#include "stdio.h"

#include "math.h"

void newton(double x[11],double y[11]); //牛顿法求解非线性方程组

void eryuanchazhi(double te[11][21],double ue[11][21]); //分片二次代数插值，求z=f(x,y) void qumianchazhi(); //曲面插值，求p(x,y)

void compare(); //比较f(x*,y*)与p(x*,y*)

void inverse(double X[10][10],int N); //求逆矩阵

double fanshu(double *p); //求向量的无穷范数

double te[11][21]={0}; //存储非线线方程组的解t

double ue[11][21]={0}; //存储非线线方程组的解u

double ze[11][21]={0}; //存储分片二次插值的解z

double x[11]={0}; //x

double y[21]={0}; //y

double inv[10][10]={0}; //存储某矩阵对应的逆矩阵

double C[10][10]={0}; //存储曲面插值的系数矩阵

double rs=0; //存储k值，用于比较函数中使用

void main()

{ int i,j;

for(i=0;i<11;i++)

for(j=0;j<21;j++)

{

x[i]=0.08*i;

y[j]=0.5+0.05*j;

}

newton(x,y);

eryuanchazhi(te,ue);

qumianchazhi();

compare();

}

void newton(double x[11],double y[11])

{

double t=0;

double u=0;

double w=0;

double v=0;

double dF[5][5]={0};

double F[5]={0};

double d[5]={0};

int i,j,k;

int ik;

int cnt,count;

double temp=0;

double m=0;

double jie[4]={0};

for(i=0;i<11;i++)

for(j=0;j<21;j++)

{ t=1;

u=1;

w=1;

v=1; //初始化

//===============求解非线性方程组=============//

do

{ dF[1][1]=-0.5*sin(t);

dF[1][2]=1;

dF[1][3]=1;

dF[1][4]=1;

dF[2][1]=1;

dF[2][2]=0.5*cos(u);

dF[2][3]=1;

dF[2][4]=1;

dF[3][1]=0.5;

dF[3][2]=1;

dF[3][3]=-sin(v);

dF[3][4]=1;

dF[4][1]=1;

dF[4][2]=0.5;

dF[4][3]=1;

dF[4][4]=cos(w);

F[1]=-(0.5*cos(t)+u+v+w-x[i]-2.67);

F[2]=-(t+0.5*sin(u)+v+w-y[j]-1.07);

F[3]=-(0.5*t+u+cos(v)+w-x[i]-3.74);

F[4]=-(t+0.5*u+v+sin(w)-y[j]-0.79);

//对Newton法的系数矩阵进行赋值,注意这里是-F() //=========================列主元高斯消元法求解方程组====================//

for (k=1;k<=3;k++)

{

ik=k;

for(cnt=k;cnt<=4;cnt++)

{

if (fabs(dF[cnt][k])>fabs(dF[ik][k]))

{

ik=cnt;

}

}