数值分析作业.docx

1、数值分析作业数值分析作业（鲁鹏：2009020979岩土工程）1已知A和b 1.1用超松弛法求解Ax=b(取松弛因子=1.4，=0，迭代9次)解：1.1.1原理：其基本思想是在高斯方法已求出x(m),x(m-1)的基础上，经过重新组合的新序列，而此新序列收敛速度加快。其算式是：xi(m)=(1-)xi(m-1)+(bijxi(m)+ xj(m-1)+gi)其中是超松弛因子，当1时，可以加快收敛速度。 1.1.2.程序：#include #include #include#define sorc 1.4#define A 9void main()int i,j,k;double aaAA=12.

2、38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743, 2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124, -1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103, 1.112336,-1.012345,3.123848,27.108437,4.10101

3、1,-3.741856,2.101023,-0.71828,-0.037585, -0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317, 0.71819,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417, 1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474, 3

4、.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782, -2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001;double bbA=2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392;double xlA,x9,tem=0.0;clrscr(

5、);for(i=0;iA;i+)xli=0;xi=0;for(k=0;k9;k+)for(i=0;iA;i+)tem=0;for(j=0;jA;j+)if(j!=i)tem=tem-aaij*xlj;xi=(tem+bbi)/aaii;xi=xli+sorc*(xi-xli);xli=xi;printf(nnttthe results are shown below!nnn);for(i=0;iA;i+)printf(ttx%d=%lfn,i+1,xi);printf(nnttpress any key to exit!n);getch();1.1.3输出结果：1.1.4讨论： 超松弛法，对G

6、-S 方法每一步得到的结果与前一步所得到的结果进行了处理，迭代收敛更快，实践证明，选择恰当的松弛系数不仅可以加快收敛速度，而且可以得到满足精度的解答！1.2 请用列主元素消去法求解Ax=b：1.2.1原理：对矩阵作恰当的调整，选取绝对值尽量大的元素作为主元素。然后进行行变换，把矩阵化为上三角阵，再进行回代，求出方程的解。1.2.2程序： #include #include #includevoid main() int i,j,k,l,m,n,r,maxi;double tem,li,x9;double aa99=12.38412,2.115237,-1.061074,1.112336,-0.

7、113584,0.718719,1.742382,3.067813,-2.031743,2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124 ,-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585

8、 , -0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317, 0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417, 1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474, 3.067813,1.112348,0.336993,-0.71828,2.

9、121314,-1.103456,3.123789,30.719334,4.446782, -2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001;double bb9=2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392;clrscr();for(i=0;i8;i+) maxi=i; for(j=i;jabs(aamaxii) for(k=i;k9;k

10、+) tem=aaik; aaik=aamaxik; aamaxik=tem; for(l=i+1;l9;l+) li=aali/aaii*(-1); for(m=i;m=0;i-)tem=0;for(j=8;ji;j-)tem=tem+aaij*xj;xi=(bbi-tem)/aaii;for(i=0;i9;i+)printf(ttx%d=%lfn,i,xi);getch();1.2.3结果输出：1.2.4讨论：列主元素消去法减小了误差传递，结果更为精确！3已知函数如下表：x12345f(x)00.693147181.09861231.36829441.6094378x678910f(x)1

11、.79175951.94591012.0794452.19722462.3025851f(x)f(1)=1f(10)=0.1试用三次样条差值求f(4.563)及f(4.563)的近似值3.1原理：根据对任意分划的三弯矩插值法：把区间分为n个点，,取区间,在这个区间上,s(x)是三次多项式，s”(x)是一次多项式，设(待定常数),则s”(x)是过（），（）两点的直线。则 （记） 将s”(x)积分两次，利用,确定任意常数，得： 记边界条件： 则由于x的划分为均匀的，则 , (其中对于来说 j=18) 则矩阵可化为：3.2程序：#include#include#include#define N 10

12、void main()int i;double x1,y1,mj1;double xN=1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,hN-1,mjN,yd0=1,ydn=0.1,dN,aN-1,bN,cN-1,ldN-1,bdN,rN-1,dtN-1;double yN=0.0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851;clrscr(); for(i=0;iN;i+) bi=2; for(i=0;iN-1;i+) hi=

13、xi+1-xi;d0=6/h0*(y1-y0)/h0-yd0); for(i=1;iN-1;i+) di=6/(hi-1+hi)*(yi+1-yi)/hi-(yi-yi-1)/hi-1);dN-1=6/hN-2*(ydn-(yN-1-yN-2)/hN-2); for(i=0;iN-1;i+) ai=hi/(hi+hi+1); ci=1-ai; bd0=b0; for(i=0;iN-1;i+) ri=ci; ldi=ai/bi; bdi+1=bi+1-ri*ldi; dt0=d0; for(i=1;i=0;i-) mji=(dti-ri*mji+1)/bdi;x1=4.563;printf(nn

14、t the results are shown below!nn); for(i=0;i=xi&x1=xi+1) y1=mji*pow(xi+1-x1),3)/(6*hi)+mji+1*pow(x1-xi),3)/(6*hi)+(yi-mji*pow(hi,2)/6)*(xi+1-x1)/hi+(yi+1-mji+1*pow(hi,2)/6)*(x1-xi)/hi; mj1=pow(x1-xi),2)/(2*hi)*mji+1-mji*pow(xi+1-x1),2)/(2*hi)+(mji-mji+1)*hi/6+(yi+1-yi)/hi; printf(t f(%lf)=%lfntdf(%l

15、f)=%lfn,x1,y1,x1,mj1); break; printf(ntplease press anykey to exit the programme!n);getch();3.3结果输出：3.4讨论：其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中不取1,x,x2,x3,(x-xj)+3为基函数,而取B样条函数3(x-xj/h)为基函数.由于三次样条函数空间是N+3维的,故我们把分点扩大到X-1,XN+1,则任意三次样条函数可用3(x-xj/h)线性组合来表示 S(x)= cj3(x-xj/h) 这样对不同插值问题,若能确定cj由解的唯一性就能求得解S(x). 样条插值因

16、为没有Runge现象，计算比高次插值函数简单，效果比高次插值函数好。且光滑性比低次分段插值函数好。它集中了高次插值函数，低次插值函数，低次分段插值函数的优点。4.用牛顿法求方程：在（0.1，1.9）中的近似根（初始近似值取为区间断点，迭代6次货误差小于0.00001）4.1原理：设函数在有限区间a,b上二阶导数存在,且满足条件） f(a)f(b)0）f”(x)在区间a,b上不变号.） f(x)0） |f(c)|/b-a|f(c)|其中c是a,b中使min|f(a),f(b)达到的一个,则对任意时近似值x0a,b,Newton迭代过程为 xk+1=(xk)=x-f(xk)/f(xk),k=1,2

17、,3用Newton法求7次方程的根，当N=1时，有： 满足Newton法的条件,所以设其初始值为1.9，故方程为：x7-28*x4+14=0 , Newton法迭代公式为： 如此一次一次的迭代，逼近x的真实根。当前后两个的差=时，就认为求出了近似的根。本程序用Newton法求代数方程(最高次数不大于10)在（a,b）区间的根4.2程序：#include #include #includedouble f1(double x)double y;y=pow(x,7)-28*pow(x,4)+14;return y;double f2(double x)double y;y=7*pow(x,6)-2

18、8*4*pow(x,3);return y;void main()int k=0;double x1,x2,f1,f2,tep, i;clrscr();printf(nntthe whole process is shown below:n);for(i=0.1;i0.00001); printf(tx0=%1.1lf-x=%lf ,i,x2); k=k+1; if(k%2=0) printf(n); getch(); 4.3结果输出：4.4讨论：由以上不同初值的迭代结果可见，取不同的初值，该迭代法会收敛于不同的解.5, 用Romberg算法求（允许误差=0.00001）5.1原理：当就停止运

19、算。5.2程序：#include #include#include#define A 1#define B 3double f(double x) double y; y=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(pow(x,2); return y; void main() int i,k=1,j; double t1010,tep; clrscr(); t00=(B-A)/2*(f(A)+f(B); t01=(t00+(B-A)*f(B+A)/2)/2; t10=(4*t01-t00)/3; printf(nntthe process is shown below!n

20、tt00=%lfn,t00); printf(tt10=%lfn,t10); while(fabs(tk0-tk-10)0.00001) k=k+1; tep=0; for(j=1;j=pow(2,k-1);j+) tep=tep+f(A+(2*j-1)*(B-A)/pow(2,k); t0k=(t0k-1+(B-A)/pow(2,k-1)*tep)/2; for(i=1;i=k;i+) tik-i=(pow(4,i)*ti-1k-i+1-ti-1k-i)/(pow(4,i)-1); printf(tt%d%d=%lfn,k,0,tk0); printf(ntpress any key to

21、exit!n);getch();5.3计算过程输出：5.4讨论：Romberge算法的优点是:1. 把积分化为代数运算,而实际上只需求T1(i),以后由递推法可得.2. 算法简单且收敛速度快,一般4或5次即能达到要求.3. 节省存储量.Romberge算法的缺点是:1. 对函数的光滑性要求较高.2. 计算新分点的数值时,这些数值的个数将成倍增加.6. 用定步长四阶Runge-Kutta法求解h=0.0005,打印(0.025), (0.045)，（0.085），（0.1），（ i=1，2，3）6.1原理： 表示 6.2程序：#include #include #include #define

22、H 0.0005double f(int n,double x,double y1,double y2,double y3)double tem;if(n=0)tem=1+0*x+0*y1+0*y2+0*y3;if(n=1)tem=0+0*x+0*y1+0*y2+y3;if(n=2)tem=1000+0*x+0*y1-1000*y2-100*y3;return tem;void main()int i,j;double yy3=0,0,0;double k34,xx=0.0005,tt3=0,0,0,tem;clrscr();printf(nnnttt the results are show

23、n below!nnn);while(1)for(i=0;i3;i+) ki0=f(i,xx,tt0,tt1,tt2);tti=yyi+H*ki0; ki1=f(i,xx+H/2,tt0,tt1,tt2);tti=yyi+H*ki1/2; ki2=f(i,xx+H/2,tt0,tt1,tt2);tti=yyi+H*ki2; ki3=f(i,xx+H/2,tt0,tt1,tt2); for(j=0;j3;j+) ttj=yyj; xx=xx+H;for(i=0;i3;i+)yyi=yyi+H/6*(ki0+2*ki1+2*ki2+k03);tti=yyi;if(fabs(xx-0.025)=0.

24、00005|fabs(xx-0.045)=0.000005|fabs(xx-0.085)=0.000005|fabs(xx-0.1)0.1)break;getch();6.3输出结果：6.4讨论：Runge_Kutta方的优点:1. 精度高,不必用别的方法求开始几点的函数值.2. 可根据f(t,y)变化的情况与需要的精度自动修改步长.3. 方法稳定.4. 程序简单,存储量少.Runge_Kutta方的缺点:每步要计算函数值f(t,y)四次,在f(t,y)，较复杂时,工作量大, 且每一步缺乏可靠的检查.另外，仅就本程序而言，在写函数的时候另外引入了一个整形变量，方便了循环体中套用，是本人的创新之处！学生：鲁鹏二零零九年十二月廿七日