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

牛人报告,供大家参考,非常鄙视抄袭者。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 …… 此处隐藏：4560字，全部文档内容请下载后查看。喜欢就下载吧 ……