数值分析编程题c语言

1、 数值分析实习报告- 12 -上机实习题一一、题目：已知A与b12.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.01

2、2345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585 A= -0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784137 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,

3、14.7138465,3.123789,-2.213474 3.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.00001b=2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.67843921用Househol

4、d变换，把A化为三对角阵B（并打印B）。2用超松弛法求解BX=b（取松弛因子=1.4,X(0)=0,迭代9次）。3用列主元素消去法求解BX=b。二、解题方法的理论依据: 1 、用Householder变换的理论依据 1 令A0=A,a(ij)1=a(ij),已知Ar_1即Ar_1=a(ij)r 2 Sr=sqrt(pow(a,2) 3 a(r)=Sr*Sr+abs(a(r+1,r)*Sr 4 y(r)=A(r_1)*u®/a® 5 Kr=(/2)*Ur的转置*Yr/a® 6 Qr=Yr-Kr*Ur 7 Ar=A(r-1)-(Qr*Ur的转置+Ur*Qr的转置) r

5、=1,2,n-22 、用超松弛法求解其基本思想：在高斯方法已求出x(m),x(m-1)的基础上，组合新的序列，从而加快收敛速度。其算式：其中是超松弛因子，当>1时，可以加快收敛速度3 、用消去法求解用追赶消去法求Bx=b的方法： , , , , q10=0 , u10=0 , x9=u19三、1计算程序：#include "math.h"#include "stdio.h"#define ge 8void main()int sign(double x);double a9=12.38412,2.115237,-1.061074,1.112336,

6、-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.101011,-3.741856,2.101023,-0.71828,-0.037

7、585,-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.713846,3.123789,-2.213474,3.067813,1.112348,0.336993,-0.71828,2.121

8、314,-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 k,h,s,w;int i,j,n,m,g;double u9,x19,y9,q9,b1910,x9;double b9=2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392; for(j=0;j&l

9、t;7;+j) /*Household 变换 */s=0.0; for(i=j+1;i<9;+i)s=s+aij*aij;s=sqrt(s);h=(aj+1j>0)?(s*s+s*aj+1j):(s*s-s*aj+1j);for(g=0;g<9;+g)if (g<=j)ug=0;else if (g=j+1) ug=aj+1j+s*sign(aj+1j);else ug=agj;for(m=0;m<9;+m)ym=0;for(n=0;n<9;+n)ym=ym+amn*un; ym=ym/h;k=0;for(i=0;i<9;+i)k=k+ui*yi;k=

10、0.5*k/h;for(i=0;i<9;+i)qi=yi-k*ui;for(n=0;n<9;+n)for(m=0;m<9;+m)amn=amn-(qm*un+um*qn);printf("Household:n");for(i=0;i<9;+i)for(j=0;j<9;+j)if (j%9=0)printf("n");printf("%-9.5f",aij);printf("n");w=1.4; /*超松弛法*/for(i=0;i<9;i+)x1i=0;for(i=0;i<

11、9;i+)for(j=0;j<9;j+)if(i=j)b1ij=0;else b1ij=-aij/aii;for(i=0;i<9;i+)b1i9=bi/aii;for(n=0;n<9;n+)for(i=0;i<9;i+)s=0;for(j=0;j<9;j+)s=s+b1ij*x1j;s=s+b1i9;x1i=x1i*(1-w)+w*s;for(i=0;i<9;i+) if (i=5)printf("n");printf("x%d=%-10.6f",i,x1i); printf("n");u0=a00

12、; /*以下是消去法*/y0=b0;for(i=1;i<9;i+)qi=aii-1/ui-1;ui=aii-qi*ai-1i;yi=bi-qi*yi-1;xge=yge/uge;for(i=ge-1;i>=0;i-)xi=(yi-aii+1*xi+1)/ui;for(i=0;i<9;i+)if (i=5)printf("n");printf(" x%d=%-10.6f",i,xi); int sign(double x) int z; z=(x>=(1e-6)?1:-1); return(z); 2打印结果：Household:1

13、2.38412 -4.89308 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000-4.89308 25.39842 6.49410 0.00000 0.00000 0.00000 0.00000 0.00000 0.000000.00000 6.49410 20.61150 8.24393 0.00000 0.00000 0.00000 0.00000 0.000000.00000 0.00000 8.24393 23.42284 -13.880070.00000 0.00000 0.00000 -0.000000.00000 0

14、.00000 0.00000 -13.8800729.69828 4.53450 0.00000 0.00000 0.000000.00000 0.00000 0.00000 0.00000 4.53450 16.00612 4.88144 0.00000 0.000000.00000 0.00000 0.00000 0.00000 0.00000 4.88144 26.01332 -4.50363 -0.000000.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -4.50363 21.25406 4.504500.00000 0.00000 0

15、.00000 -0.00000 0.00000 0.00000 -0.00000 4.50450 14.53412x0=1.073409 x1=2.272579 x2= -2.856601x3=2.292514 x4=2.112165 x5= -6.422586x6=1.357802 x7=0.634259 x8= -0.587042x0=1.075799 x1=2.275744 x2= -2.855515x3=2.293099 x4=2.112634 x5= -6.423838x6=1.357923 x7=0.634244 x8= -0.587266四、问题讨论：此程序具有很好的通用性。在G

16、S方法的基础上，已经求出x的第m解，第m-1基础上，经过重新组合得新的序列，而在此新序列收敛速度加快。上机实习题二一、题目：已知函数值如下表：x12345f(x)00.69314781.09861231.38629441.6094378x678910f(x)1.7917951.94591012.0794452.19722462.3025851f(x)f(1)=1f(0)=0.1试用三次样条插值求f(4.563)及f(4.563)的近似值。二、解题方法的理论依据：任意划分的三弯矩插值法以及方程组解法中的三对角阵追赶算法。应用三次样条插值法能够对函数产生很好的逼近效果。而追赶算法又具有计算量少、方

17、法简单、算法稳定的特点。方法应用条件：适用于求复杂函数在给定区间内某一点的函数值，给出函数f(x)在区间a,b中的n个插值点，并且给出函数在区间端点处的值。三、1计算程序：#include "stdio.h"#include "math.h"#define n 11#define ge 11void main()int i,m;double e,s,E,q12,u12,y12,c12,w12;double b12=2,0,4.15888308,6.5916738,8.3177664,9.6566268,10.750557,11.6754606,12.47

18、667,13.1833476,13.8155106,14.0155106;double a1212=-2,-4,0,0,0,0,0,0,0,0,0,0,0;ann-1=4;ann=2;for(i=1;i<11;i+)aii-1=1;aii=4;aii+1=1;u0=a00; /*消去法求ci*/y0=b0;for(i=1;i<12;i+)qi=aii-1/ui-1;ui=aii-qi*ai-1i;yi=bi-qi*yi-1;cge=yge/uge;for(i=ge-1;i>=0;i-)ci=(yi-aii+1*ci+1)/ui; printf("请输入要插的值:&q

19、uot;);scanf("%lf",&E); for(i=0;i<12;i+) e=fabs(E-i); if(e>=2) wi=0;else if(e<=1) wi=0.5*fabs(e*e*e)-e*e+2.0/3.0;elsewi=(-1.0/6.0)*fabs(e*e*e)+e*e-2*fabs(e)+4.0/3.0; s=0.0;for(i=0;i<12;i+)s=s+ci*wi;printf("f(%lf)=%lf",E,s);printf("n");printf("请输入要求的导

20、数的值:");scanf("%d",&m);printf("f(%d)=%lfn",m,(cm+1-cm-1)/2.0); 输出结果：请输入要插的值: 4.563f(4.563)=1.517932请输入要求的导数的值: 4.563f(4.563)= 0.249350四、问题讨论：在给均匀分划的插值函数x赋值时，由于使用for循环，误将xi=i+1写成xi=i，导致运算错误。此程序具有一定通用性，对于任意划分的三弯矩插值法，只许改动xi即可。求解方程组Mj时，要用到三对角方程组的追赶法（也称Thomas算法）。变量较多，应注意区分。求导时

21、注意正负号。上机实习题三一、题目： 用Newton法求方程 X7-28x4+14=0在（0.1,1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。二、解题方法及理论依据：Newton迭代法是平方收敛于方程f(x)=0在区间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,3算

22、法：令故以1.9为起点三、1计算程序：#include "math.h"main() float x,y,f,f1; scanf("%f",&x); do y=x; f=pow(y,7)-28*pow(y,4)+14; /*定义f(x)的表达式*/ f1=7*pow(y,6)-112*pow(y,3); /*定义f(x)的表达式*/ x=y-f/f1; /* Newton迭代法*/ while(fabs(x-y)>=1e-5); /*控制误差小于0.00001*/ printf("nThe result of the questi

23、on is %fn",x);2打印结果：请输入端点值：1.9 0.1x=0.845497 3.030577四、问题讨论：程序较为简单。它的几何意义为xk+1是f(x)在点xk的切线与x轴交点，故也称为切线法，它是平方收敛的，此处取xk=1.9收敛性较好，要注意判断f(xk)是否为零。上机实习题四一、题目：用Romberg算法求 (允许误差=0.00001)。二、解题方法及理论依据：龙贝格（Romberg）方法求数值积分T1(0)=(b-a)/2*f(a)+f(b)T1(l)=（1/2）* T1(l-1)+(b-a)/2l-1*fa+(2i-1)*(b-a)/2lTm+1k-1=4mT

24、m(k)-Tm(k-1)/(4m-1)三、1计算程序：#include "math.h"int a=1,b=3;double f(double x) /*求f(x)=3xx1.4(5x+7)sinx2的值*/ double z; z=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return (z);double s(int l) /*求T1(l) 中的(b-a)/2l-1*fa+(2i-1)*(b-a)/2l*/ extern a,b; int i,m; double z=0.0; m=pow(2,l-1); for(i=1;i<=m;

25、i+)z+=f(a+1.0*(2*i-1)/m); z*=1.0*(b-a)/m; return (z);main() extern a,b; double T2020; int m,n,l=0; T10=(b-a)/2.0*(f(a)+f(b); do /*龙贝格（Romberg）算法*/ l+; T1l=0.5*(T1l-1+s(l); n=l-1; for(m=2;n>=0;m+,n-) /*求解Tl(0)*/ Tmn=(pow(4,m-1)*Tm-1n+1-Tm-1n)/(pow(4,m-1)-1.0); while(fabs(Tl0-Tl-10)>=1e-5); prin

26、tf("nT%d0=%f",l,Tl0);2打印结果：T80=440.536017四、问题讨论：此程序较繁，计算T1(k)需要复化梯形公式，还要用到Richardson外推法，构造新序列，计算新分点的值时，这些数值个数成倍增加。应用给出所要求的误差，当|Tl(0)-Tl+1(0)|<时控制循环。程序具有广泛的通用性。上机实习题五一、题目：用定步长四阶Runge-Kutta法求解 dy1/dt=1 dy2/dt=y3 dy3/dt=1000-1000y2-100y3 y1(0)=0 y2(0)=0 y3(0)=0h=0.0005,打印yi(0.025), yi(0.04

27、5), yi(0.085), yi(0.1),(i=1,2,3)二、解题方法及理论依据：高阶方程组的Runge-Kutta解法 Yn+1=Yn+(1/6)*(K1+2K2+2K3+K4) K1=h*F(xn,Yn) K2=h*F(xn+h/2,Yn+K1/2) K3=h*F(xn+h/2,Yn+K2/2) K4=h*F(xn+h,Yn+K3) 适用条件：使用于那些用普通的积分方法解不了的微分方程组.只要知道函数之间的关系和初值就可以不用解出表达式而直接求解函数在要求点的值。三、1计算程序：#include <stdio.h>main() double t,h=0.0005,y1=0

28、,y2=0,y3=0,ky15,ky25,ky35; for(t=h;t<=0.1001;t+=h) /*Runge-Kutta算法具体过程*/ ky11=h*1; ky21=h*y3; ky31=h*(1000-1000*y2-100*y3); ky12=h*1; ky22=h*(y3+ky31/2); ky32=h*(1000-1000*(y2+ky21/2)-100*(y3+ky31/2); ky13=h*1; ky23=h*(y3+ky32/2); ky33=h*(1000-1000*(y2+ky22/2)-100*(y3+ky32/2); ky14=h*1; ky24=h*(y3+ky33); ky34=h*(1000-1000*(y2+ky23)-100*(y3+