《北航数值分析第三次大作业》由会员分享，可在线阅读，更多相关《北航数值分析第三次大作业（25页珍藏版）》请在金锄头文库上搜索。
1、数值分析第三次大作业姓名:袁二凯学号:SY0917145联系方式:2011年 12月 8日一、 算法的设计方案：（一） 、总体方案设计：（1）解非线性方程组。将给定的 当作已知量代入题目给定的非线性方程组，(,)ixy求得与 相对应的数组 tij,uij。(,)ixy（2）分片二次代数插值。通过分片二次代数插值运算，得到与数组 t1121,u1121对应的数组 z1121，得到二元函数 z= 。(,)ifxy（3）曲面拟合。利用 xi,yj,z1121建立二维函数表，再根据精度的要求选择适当 k值，并得到曲面拟合的系数矩阵 Crs。（4）观察和 的逼近效果。观察逼近效果只需要重复上面（1）和（
2、2）的过(,)ipxy程，得到与新的插值节点 对应的 ，再与对应的 比较即可，这里i (,)ifxy(,)ipxy求解 可以直接使用（3）中的 Crs和 k。(,)ixy（二）具体算法设计：（1）解非线性方程组牛顿法解方程组 的解 ，可采用如下算法：()0Fx*1）在 附近选取 ，给定精度水平 和最大迭代次数 M。*x()D02）对于 执行,1kML 计算 和 。()kFx()k 求解关于 的线性方程组()()()()kkxFx 若 ，则取 ，并停止计算；否则转。()()k*() 计算 。(1)()()kkkxx 若 ，则继续，否则，输出 M 次迭代不成功的信息，并停止计算。M（2）分片双二次
3、插值给定已知数表以及需要插值的节点，进行分片二次插值的算法：设已知数表中的点为： ，需要插值的节点为0(,1),ijxihnyjmL。(,)xy1) 根据 选择插值节点 ：(,)xy(,)ijxy若 或 ，插值节点对应取 或 ，12h12nh1in若 或 ，插值节点对应取 或 。yyjim若,2,2iijjhxxinyyjm则选择 为插值节点。(,)1,;1,)krxirj2）计算11()(,)()(1,)itktikjtrtjrxliyl j%插值多项式的公式为：1(,)()(,)jikrkrkirjpxylxyf注：本步进行插值运算的是 ，利用 与 的对应关系就可以得到,tu,ij,tu与
4、 的对应关系。z(,)ijxy（3）曲面拟合根据插值得到的数表 进行曲面拟合的过程：,(,)ijijxyf1） 根据拟合节点和基底函数写出矩阵 B 和 G：0101101()()()()kknnxxBLMO0101101()()()()kkmmyyLMO2） 计算 。1TTCBUG在这里，为了简化计算和编程、避免矩阵求逆，记：，1()TA1()TTD对上面两式进行变形，得到如下两个线性方程组：，()TTBAU()TTGD通过解上述两个线性方程组，则有： TCA3） 对于每一个 ， 。(,)ijxy*0(,)()krsij sijrpxyxy4） 拟合需要达到的精度条件为：。*270(,)1nm
5、ijijijxyu其中 对应着插值得到的数表 中 的值。iju,(,)ijijf(,)ijfxy5） 让 k 逐步增加，每一次重复执行以上几步，直到成立。此时的 k 值就是要求解最小的 k。*270(,)10nmijijijpxyu二、 源程序：#include#include#include #include #include#include #define Epsilon1 1e-12 /*解线性方程组时近似解向量的精度*/#define M 200 /*解线性方程组时的最大迭代次数*/#define N 10 /*求解迭代次数时假设的 k的最大值，用于定义包含 k的存储空间*/void
6、Newton(); /*牛顿法求解非线性方程组子程序*/void fpeccz(); /*分片二次代数插值子程序*/void qmnh(); /*曲面拟合子程序*/void duibi(); /*对比和 p逼近效果的子程序*/double x11,y21,t1121,u1121;/*定义全局变量*/double z1121,C1010;double kz;void Newton(double x11,double y21)/*牛顿法求解非线性方程组子程序*/ double X4,dx4,F4,dF44,temp,m,fx,fX;int i,j,k,l,p,ik,n; for(i=0;i=0;k
7、-)temp=0;for(l=k+1;l0.7) sij=4;else for(m=2;m0.2*m-0.1)&(tij1.4) rij=4; elsefor(m=2;m0.4*m-0.2)&(uij1e-7) for(i=0;i=0;m-)temp=0;for(i=m+1;i=0;m-)temp=0;for(i=m+1;i=k;i+) temp+=Dil*GtG1mi; Dml=(Gtml-temp)/GtG1mm;for(i=0;i=k;i+)for(j=0;j=k;j+)temp=0; for(m=0;m=20;m+) temp+=Aim*Djm; Cij=temp;sigma=0;/*
8、归零，开始计算 sigma值*/for(i=0;i=10;i+)for(j=0;j=20;j+)pij=0;for(i1=0;i1=k;i1+)for(j1=0;j1=k;j1+) pij+=Ci1j1*pow(xi,i1)*pow(yj,j1); sigma+=pow(pij-zij,2);printf(k=%d sigma=%.12en,k,sigma); k=k+1;k-;printf(达到精度要求时的 k和 sigma值：n);printf(k=%d sigma=%.12en,k,sigma); kz=k; /*定义为全局变量，便于 duibi()调用*/void duibi() /*
9、对比和 p逼近效果的子程序*/ int i,j,i1,j1;double p85;for(i=0;i=7;i+) for(j=0;j=4;j+)xi=0.1*(i+1);yj=0.5+0.2*(j+1); /*重新输入节点*/Newton(x,y);fpeccz(t,u); /*求解（ x*，y *）*/ for(i=0;i=7;i+) /*求解 p（x *，y *）*/for(j=0;j=4;j+)pij=0;for(i1=0;i1=kz;i1+) for(j1=0;j1=kz;j1+) pij+=Ci1j1*pow(xi,i1)*pow(yj,j1); printf(x%d=%.6f y%
10、d=%.6fn,i+1,xi,j+1,yj);printf(f(x%d,y%d)=%.12en,i+1,j+1,zij);printf(p(x%d,y%d)=%.12en,i+1,j+1,pij);printf(=%.12en,zij-pij); /*数表 xi*,yi*, (xi*,yi*),p（x i*,yi*）*/void main() int i,j;for(i=0;i=10;i+) for(j=0;j=20;j+)xi=0.08*i;yj=0.5+0.05*j; /*输入节点*/Newton(x,y);fpeccz(t,u);for(i=0;i=10;i+) for(j=0;j=20
11、;j+)printf(x%d=%.6f y%d=%.6f z%d%d=%.12e n,i,xi,j,yj,i,j,zij); /*数表：x i,yi, (xi,yi)*/qmnh();for(i=0;i=kz;i+) for(j=0;j=kz;j+) printf(C%d%d=%.12e n,i,j,Cij);/*数表：Crs*/duibi();三、运行结果输出1、数表数表：x i,yi,(xi,yi)x0=0.000000 y0=0.500000 z00=4.465040184807e-001x0=0.000000 y1=0.550000 z01=3.246832629277e-001x0=
12、0.000000 y2=0.600000 z02=2.101596866827e-001x0=0.000000 y3=0.650000 z03=1.030436083160e-001x0=0.000000 y4=0.700000 z04=3.401895562675e-003x0=0.000000 y5=0.750000 z05=-8.873581363800e-002x0=0.000000 y6=0.800000 z06=-1.733716327497e-001x0=0.000000 y7=0.850000 z07=-2.505346114666e-001x0=0.000000 y8=0.9
13、00000 z08=-3.202765063876e-001x0=0.000000 y9=0.950000 z09=-3.826680661097e-001x0=0.000000 y10=1.000000 z010=-4.377957667384e-001x0=0.000000 y11=1.050000 z011=-4.857589414438e-001x0=0.000000 y12=1.100000 z012=-5.266672548835e-001x0=0.000000 y13=1.150000 z013=-5.606384797965e-001x0=0.000000 y14=1.200000 z014=-5.877965387677e-001x0=0.000000 y15=1.250000 z015=-6.082697790899e-001x0=0.000000 y16=1.300000 z016=-6.221894528764e-001x0=0.000000 y17=1.350000 z017=-6.296883781856e-001x0=0.000000 y18=1.400000 z018=-6.308997600028e-001x0=0.000000 y19=1.4
