《研究生数值分析课后题(上机编程).doc》由会员分享，可在线阅读，更多相关《研究生数值分析课后题(上机编程).doc（14页珍藏版）》请在三一办公上搜索。
1、 2009级研究生数值分析上机作业院系 电气工程学院专业 控制理论与控制工程姓名 马凯学号 2009020954指导教师 代新敏2009年12月29日第一题（二问）：超松弛法求方程组根1.解题理论依据或方法应用条件:超松弛算法是在GS方法已求出x(m),x(m-1)的基础上，经过重新组合得到新序列。如能恰当选择松弛因子，收敛速度会比较快。当1时，称为超松弛法，可以用来加速收敛。其具体算法为： 2.计算程序（使用软件：VC）：#include#define w 1.4main()float a1010= 0,0,0,0,0,0,0,0,0,0,0,12.38412,2.115237,-1.061
2、074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743,0,2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124,0,-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,0,1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.1010
3、23,-0.71828,-0.037585,0,-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317,0,0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417,0,1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474,0,3.067813,1.112
4、348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782,0,-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001; float b101=0,2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392; float x1010=0,0,0,0,0,0,0,0,0,0; /*由x
5、(0)=0得到其第一列全为零*/ float sum1=0,sum2=0; int i,m,j; for(m=1;m=9;m+) for(i=1;i=9;i+) sum1=0; for(j=1;j=(i-1);j+)sum1+=(-aij/aii)*xjm; /*计算第一个累加和*/ sum2=0; for(j=(i+1);j=9;j+)sum2+=(-aij/aii)*xjm-1; /*计算第二个累加和*/ xim=(1-w)*xim-1+w*(sum1+sum2+bi0/aii); /*用SOR方法计算*/ printf(x1为：%lfn,x19); printf(x2为：%lfn,x29
6、); printf(x3为：%lfn,x39); printf(x4为：%lfn,x49); printf(x5为：%lfn,x59); printf(x6为：%lfn,x69); printf(x7为：%lfn,x79); printf(x8为：%lfn,x89); printf(x9为：%lfn,x99);3.计算结果4.问题讨论（误差分析、上机出现情况等）这道题目是所有题目中编写最顺利的，一次即顺利得出结果，当然这道题目还是有应该注意到地方，一是注意两个求和的清零，二是注意下标，不要弄混行标和列标。第一题（三问）：列主元素法求方程组根1.解题理论依据或方法应用条件: 所谓列主元消去法是，
7、对矩阵作恰当的调整，选取绝对值最大的元素作为主元素。然后把矩阵化为上三角阵，再进行回代，求出方程的解。算法为： , , , , q10=0 , u10=0 , x9=u192.计算程序（使用软件：VC）：#include#includemain()double a1011= 0,0,0,0,0,0,0,0,0,0,0,0,12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743,2.1874369,0,2.115237,19.141823,-3.125432,-1.012345,2.1
8、89736,1.563849,-0.784165,1.112348,3.123124,33.992318,0,-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,-25.173417,0,1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585,0.84671695,0,-0.113584,2.189736,2.031454,4.101011,19.897918,0
9、.431637,-3.111223,2.121314,1.784317,1.784317,0,0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417,-86.612343,0,1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474,1.1101230,0,3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3
10、.123789,30.719334,4.446782,4.719345,0,-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001,-5.6784392; double x10,y11; double sum=0,max; int i,j,j1,m,n,k,g,h; for(j=1;j=8;j+)max=ajj; for(k=j+1;kfabs(max)max=akj; for(k=j;k=9;k+) if(akj/max=1)g=k; /*确定本列绝对值最大元素所在行：g行*/
11、for(h=j;h=10;h+) yh=ajh; for(h=j;h=10;h+) ajh=agh; for(h=j;h=10;h+) agh=yh; /*将g行与J行所有元素进行交换*/ for(i=j+1;i=9;i+) for(j1=j+1;j1=1;m-) sum=0; for(n=m+1;n=9;n+) sum+=xn*amn; printf(x%d为：%lfn,m,xm=(am10-sum)/amm); /*从x9到x1逐个解出结果*/3.计算结果4.问题讨论（误差分析、上机出现情况等）（1）这道题目出现的最大问题在 for(j1=j+1;j1=10;j1+)aij1-=(aij/
12、ajj)*ajj1; ，刚开始我错写为：for(j1=j;j1=10;j1+)aij1-=(aij/ajj)*ajj1; ，也就是j1应该从j+1开始取值，如果从j开始取，会导致aij1-=(aij/ajj)*ajj1; 一句中aij是零，这样后面的aij1全部都将是原来的数值，不会变化，导致错误。（2）对于列主元素消元法，个人感觉将a,b两个矩阵，放在一个矩阵中进行编程，较方便。第三题：三次样条插值求近似值1.解题理论依据或方法应用条件:依据三弯矩插值法列出相应的方程，然后将x的值带入（在求解三弯矩方程的参数时，应用追赶法求解方程组求出相应参数值）。 ，由，得出：2.计算程序（使用软件：VC
13、）：#include#includemain() int i; double s,z; /*定义函数近似值，导数近似值*/ double b11=0,2,2,2,2,2,2,2,2,2,2,p11,g11,l11; /*追赶法中变量b,l*/ double m11=0,0,0,0,0,0,0,0,0,0,0;doublef11=0,0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851; /*原插值点对应函数值*/ double y11=0,0.5,0.5,0.5,0
14、.5,0.5,0.5,0.5,0.5,0.5,1; /*变量*/ double x10=0,1,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5; /*变量*/doubled11=0,6*(f2-1),3*(f3+f1-2*f2),3*(f4+f2-2*f3),3*(f5+f3-2*f4),3*(f6+f4-2*f5),3*(f7+f5-2*f6),3*(f8+f6-2*f7),3*(f9+f7-2*f8),3*(f10+f8-2*f9),6*(0.1-f10+f9); p1=2;g1=d1; for(i=1;i=1;i-) mi=(gi-xi*mi+1)/pi; /*应用追赶法
15、求出最终方程解*/s=m4*pow(5-4.563,3)/6+m5*pow(4.563-4,3)/6+(f4-m4/6)*(5-4.563)+(f5-m5/6)*(4.563-4);z=-m4/2*pow(5-4.563,2)+m5/2*pow(4.563-4,2)-(f4-m4/6)+(f5-m5/6);printf(函数f(4.563)近似值 ：%.16lfn导数f(4.563)近似值：%.16lfn,s,z); /*求出函数值和导数值并输出*/ 3.计算结果4.问题讨论（误差分析、上机出现情况等）（1）在做数组相关的编程时，特别注意下标的使用，因为数组下标是从0开始的，以后可以统一将下标
16、为0的数组位置置0，从下标为1的位置开始赋值使用。（2）编程的过程出现问题如下：第二个for循环i-误写成i+，造成无限循环错误在乘法的运算中漏掉了两个小括号间的“”号第四题：牛顿法求方程近似根1.解题理论依据或方法应用条件:依据牛顿迭代法进行计算，设函数在有限区间a,b上二阶导数存在，且满足条件. ;. 在区间a,b上不变号；. ;. ，其中c是a,b中使达到的一个.则对任意初始近似值， Newton迭代过程为，k=0，1，2所生成的迭代序列平方收敛于方程在区间a,b上的唯一解.对本题：所以：2.计算程序（使用软件：VC）：#includevoid main()double x0,x1=1.
17、9,f0,f1;for(x0=1.0;fabs(x0-x1)=0.00001;)x0=x1; f0=pow(x0,7)-28*pow(x0,4)+14; f1=7*pow(x0,6)-112*pow(x0,3); if(f1=0)break; else x1=x0-f0/f1;if(f1=0)printf(出错n);else printf(方程的解为：%.16lfn,x1);3.计算结果4.问题讨论（误差分析、上机出现情况等）（1）在做除法运算的时候，应对除数是否等于零进行判断。（2）编程的过程出现问题如下：在编写程序过程中，不小心将输入换成了全角输入造成错误忘记头文件中增加#include，
18、导致pow函数无法使用。第五题：用Romberg算法求积分近似值1.解题理论依据或方法应用条件:依据Romberg法列出相应的方程，先求出第一列的相应元素，然后利用Romberg中三角法求出后面列的元素，将得到的第一行元素取差的绝对值，当小于误差要求的时候可以停止运算，公式为： m= 1,2,l. k=1,2, ,l-m+1 具体步骤为：（1） 先求出按梯形公式所得的积分值：（2） 把区间2等分，求出两个小梯形面积之和，记为，即，由外推法得，（3） 把区间再等分，由Richardson外推法，构造新序列，m= 1,2,l. k=1,2, ,l-m+1（4） 或就停止计算，否则回到(3),计算2
19、.计算程序（使用软件：VC）：方法一：#include#include#define f(x) pow(3,(x)*pow(x),1.4)*(5*(x)+7)*sin(x)*(x)void main()double t99; double sum=0; int j=1,i=1,m=0,k=1; double a=1,b=3;t01=(b-a)/2*(f(a)+f(b);sum=f(a+(2*1-1)*(b-a)/pow(2,1); t11=0.5*(t01+sum*(b-a)/pow(2,0); t02=(pow(4,1)*t11-t01)/(pow(4,1)-1); /*先求出t01，t02
20、的值*/ for(j=1;fabs(t0j+1-t0j)0.00001;) /*判断误差是否符合条件*/ j+; sum=0; /*累加和清零*/ for(i=1;i=1;k-) m+;tk-1m+1=(pow(4,m)*tkm-tk-1m)/(pow(4,m)-1); /*依次用三角形法求出其他元素的值*/ m=0; printf(积分的值为：%.16lfn,t0j+1);方法二：（逐列求值，较浪费电脑资源、时间）#include#include#define f(x) pow(3,(x)*pow(x),1.4)*(5*(x)+7)*sin(x)*(x)void main()double t
21、99; double sum; int j=1,i=1,m=1,k=1,l=1; double a=1,b=3;t01=(b-a)/2*(f(a)+f(b); for(l=1;l=7;l+) sum=0; for(i=1;i=pow(2,l-1);i+) sum+=f(a+(2*i-1)*(b-a)/pow(2,l); tl1=0.5*(tl-11+sum*(b-a)/pow(2,l-1); for(m=1;m=7;m+) for(k=1;k0.00001;j+); printf(积分的值为：%.16lfn,t0j+1); 3.计算结果4.问题讨论（误差分析、上机出现情况等）（1）这道题目是前
22、三道题目中耗时最长、但也是学到最多东西的题目，一开始使用的是第二个方法进行编程，VC报错，一直未找到哪里出错，最后发现原来是将int j=1,i=1,m=1,k=1,l=1; 放在了 t01=(b-a)/2*(f(a)+f(b); 之后，VC规定必须将所有的定义放在程序开头，这是自己以前所不知道的。（2）改掉上面的错误后，系统不再报错，但是数值不对，逐行看程序发现原来是sum未清零，导致后面的数据出错。（3）改掉上面的错误后，得到的数据依然不对，纠缠很久以及翻阅C语言书籍发现原来是#define f(x) pow(3,(x)*pow(x),1.4)*(5*(x)+7)*sin(x)*(x) 语
23、句中，X必须加上括号，这样才能保证运算顺序正确。（4）这样程序的结果正确了，但是第二种方法较浪费电脑资源、时间，故对原程序更改有了第一种利用循环、用多少数据求多少数据的方法。第六题：定步长四阶Runge-Kutta法求方程组1.解题理论依据或方法应用条件:依据定步长四阶Runge-Kutta法进行计算，公式为：2.计算程序（使用软件：VC）：#includemain()float y4205=0,0,0,0; /*定义第一列值全为零*/ int i,j; float k1,k2,k3,k4; float h=0.0005; for(j=1;(j*0.0005)=0.1;j+) for(i=1;
24、i=3;i+) /*本解法未使用0行，故从1开始取值*/if(i=1) k1=1;k2=1;k3=1;k4=1; y1j=y1j-1+h/6*(k1+2*k2+2*k3+k4); /*计算第一行数值*/ if(i=2) k1=y3j-1; k2=y3j-1+0.5*h*(1000-1000*y2j-1-100*y3j-1);k3=y3j-1+0.5*h*(1000-1000*(y2j-1+0.5*h*y3j-1)-100*(y3j-1+0.5*h*(1000-1000*y2j-1-100*y3j-1);k4=y3j-1+h*(1000-1000*(y2j-1+0.5*h*(y3j-1+0.5*
25、h*(1000-1000*y2j-1-100*y3j-1)-100*(y3j-1+0.5*h*(1000-1000*(y2j-1+0.5*h*y3j-1)-100*(y3j-1+0.5*h*(1000-1000*y2j-1-100*y3j-1); y2j=y2j-1+h/6*(k1+2*k2+2*k3+k4); /*计算第二行数值*/ if(i=3) k1=1000-1000*y2j-1-100*y3j-1; k2=1000-1000*(y2j-1+0.5*h*y3j-1)-100*(y3j-1+0.5*h*(1000-1000*y2j-1-100*y3j-1);k3=1000-1000*(y
26、2j-1+0.5*h*(y3j-1+0.5*h*(1000-1000*y2j-1-100*y3j-1)-100*(y3j-1+0.5*h*(1000-1000*(y2j-1+0.5*h*y3j-1)-100*(y3j-1+0.5*h*(1000-1000*y2j-1-100*y3j-1);k4=1000-1000*(y2j-1+h*(y3j-1+0.5*h*(1000-1000*(y2j-1+0.5*h*y3j-1)-100*(y3j-1+0.5*h*(1000-1000*y2j-1-100*y3j-1)-100*(y3j-1+h*k3);y3j=y3j-1+h/6*(k1+2*k2+2*k3
27、+k4); /*计算第三行数值*/ printf(y1(0.025),y2(0.025),y3(0.025)值分别为：%f,%f,%fn,y150,y250,y350); printf(y1(0.045),y2(0.045),y3(0.045)值分别为：%f,%f,%fn,y190,y290,y390); printf(y1(0.085),y2(0.085),y3(0.085)值分别为：%f,%f,%fn,y1170,y2170,y3170); printf(y1(0.100),y2(0.100),y3(0.100)值分别为：%f,%f,%fn,y1200,y2200,y3200); 3.计算结果4.问题讨论（误差分析、上机出现情况等）（1）本题的编写过程相对前一题较为顺利，过程中只出现一个较大问题，原程序中我在for循环体内使用的是if() else if() else () 的结构，无错，但运行不出结果，改用if () if () if () 后出现结果，至今仍未找清楚原因。（2）另一个小问题是，在计算第三行数据时，由于要使用前一列二行的数据，刚开始直接将复制进来，稍后发现这时的K值是二行中的数据，放在第三行中将出错，由于发现及时，未对后面程序编写造成影响。（故后续将K值都由具体的算式替代，避免出错）
