浙江大学数值分析C语言编程习题_第1页
浙江大学数值分析C语言编程习题_第2页
浙江大学数值分析C语言编程习题_第3页
浙江大学数值分析C语言编程习题_第4页
浙江大学数值分析C语言编程习题_第5页
已阅读5页,还剩39页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

1、精选优质文档-倾情为你奉上C语言编程习题第二章习题2-25. 用二分法编程求 6x4 -40x2+9=0 的所有实根。#include <stdio.h>#include <math.h>#define N 10000double A,B,C;double f(double x) return (A*x*x*x*x+B*x*x+C);void BM(double a,double b,double eps1,double eps2) int k; double x,xe;double valuea = f(a);double valueb = f(b);if (valu

2、ea > 0 && valueb > 0 | valuea <0 && valueb < 0) return;printf("Finding root in the range: %.3lf, %.3lfn", a, b); for(k=1;k<=N;k+) x=(a+b)/2; xe=(b-a)/2; if(fabs(xe)<eps2 | fabs(f(x)<eps1) printf("The x value is:%gn",x); printf("f(x)=%gnn&

3、quot;,f(x); return; if(f(a)*f(x)<0) b=x; else a=x; printf("No convergence!n");int main() double a,b,eps1,eps2,step,start; printf("Please input A,B,C:n"); scanf("%lf %lf %lf",&A,&B,&C); printf("Please input a,b, step, eps1,eps2:n"); scanf("%

4、lf %lf %lf %lf %lf",&a,&b,&step,&eps1,&eps2); for (start=a; (start+step) <= b; start += step) double left = start;double right = start + step;BM(left, right, eps1, eps2); return 0;运行:Please input A,B,C:6 -40 9Please input a,b, step, eps1,eps2:-10 10 1 1e-5 1e-5Finding roo

5、t in the range: -3.000, -2.000The x value is:-2.53643f(x)=-0.Finding root in the range: -1.000, 0.000The x value is:-0.f(x)=0.Finding root in the range: 0.000, 1.000The x value is:0.f(x)=0.Finding root in the range: 2.000, 3.000The x value is:2.53643f(x)=-0.有时若把判别语句if(fabs(xe)<eps2 | fabs(f(x)<

6、;eps1)改为if(fabs(xe)<eps2 && fabs(f(x)<eps1)会提高精度,对同一题运行结果:Finding root in the range: -3.000, -2.000The x value is:-2.53644f(x)=-4.26496e-007Finding root in the range: -1.000, 0.000The x value is:-0.f(x)=-7.3797e-006Finding root in the range: 0.000, 1.000The x value is:0.f(x)=-7.3797e-0

7、06Finding root in the range: 2.000, 3.000The x value is:2.53644f(x)=-4.26496e-007习题2-35. 请用埃特金方法编程求出x=tgx在4.5(弧度)附近的根。#include <stdio.h>#include <math.h>#define N 100#define PI 3.void SM(double x0,double eps) int k; double x; double x1,x2; for(k=1;k<=N;k+) x1=sin(x0)/cos(x0); x2=sin(x

8、1)/cos(x1); x=x0-(x1-x0)*(x1-x0)/(x2-2*x1+x0); if(fabs(x-x0)<eps) printf("The x value is:%gn",x); return; x0=x; printf("No convegence!n");int main() double eps,x0; printf("Please input eps,x0:n"); scanf("%lf %lf",&eps,&x0); SM(x0,eps); return 0;运行:P

9、lease input eps,x0:1e-5 4.5The x value is:4.49341习题2-411请编出用牛顿法求复根的程序,并求出 P(z)=z4-3z3+20z2+44z+54=0 接近于z0=2.5+4.5i 的零点。#include <cstdio>#include <cmath>#define MAX_TIMES 1000typedef struct double real, image; COMPLEX;COMPLEX Aa5=54,0,44,0,20,0,-3,0,1,0;COMPLEX Bb4=44,0,40,0,-9,0,4,0;COMP

10、LEX zero = 0, 0;double eps1=1e-6;double eps2=1e-6;COMPLEX multi(COMPLEX a,COMPLEX b) COMPLEX result; result.real = a.real * b.real - a.image * b.image; result.image = a.image * b.real + a.real * b.image; return result;COMPLEX Div(COMPLEX a,COMPLEX b) COMPLEX z3; double s; s=(b.real*b.real)+(b.image*

11、b.image); z3.real=b.real; z3.image=-b.image; z3=multi(a,z3); z3.real=z3.real/s; z3.image=z3.image/s; return z3;COMPLEX add(COMPLEX a,COMPLEX b) COMPLEX result; result.real = a.real + b.real; result.image = a.image + b.image; return result;COMPLEX subtract(COMPLEX a, COMPLEX b) COMPLEX result; result

12、.real = a.real - b.real; result.image = a.image - b.image; return result;COMPLEX times(COMPLEX z,int n) int i; COMPLEX result=1, 0; for (i=0;i<n;i+) result=multi(result,z); return result;double distance(COMPLEX a, COMPLEX b) return sqrt(a.real - b.real) * (a.real - b.real) + (a.image - b.image) *

13、 (a.image - b.image);double complex_abs(COMPLEX a) return sqrt(a.real * a.real + a.image * a.image);COMPLEX f(COMPLEX x) int i; COMPLEX result=zero; for (i=0;i<5;i+) result=add(result,multi(Aai,times(x,i); return result;COMPLEX ff(COMPLEX x) int i; COMPLEX result=zero; for (i=0;i<4;i+) result=

14、add(result,multi(Bbi,times(x,i); return result;int main() COMPLEX z0,z1,result; double x,y; int k; printf("please input x,yn"); scanf("%lf %lf",&x,&y); z0.real=x; z0.image=y; for(k=0; k<MAX_TIMES; k+) z1 = subtract(z0, Div(f(z0), ff(z0); result = f(z1); if (distance(z0

15、,z1)<eps1 | complex_abs(result)<eps2) printf("The root is: z=(%.3f + %.3fi), f(z)=(%e + %ei)n", z1.real,z1.image, result.real, result.image); return 0; z0 = z1; printf("No convergence!n"); return 0;运行:please input x,y2.5 4.5The root is: z=(2.471 + 4.641i), f(z)=(-1.e-007 +

16、-1.e-007i)习题2-62. 请编程用劈因子法求高次方程 x4+ x 3+5 x 2+4 x +4=0的所有复根。#include<cstdio>#include<cmath>int main()float b20,c20;float delta;float eps;int print=0;float fru0,fru1,s0,s1,frv0,frv1;float deltau,deltav;float u,v;float a20;int n,j;float r0,r1;float r;int i, k;printf("Please input max

17、exponent:n");scanf("%d",&n);printf("Please input epslion:n");scanf("%f",&eps);printf("Please input the coefficient:n");for(i=0;i<=n;i+) scanf("%f,",&ai);for(j=0;j<=n;j+) printf("+%.3fX%d",aj,n-j);printf("n")

18、; for(k=1;k<=2;k+) u=an-1/an-2;v=an/an-2;if(k=2) u=0;v=4;while(1) b0=a0;b1=a1-u*b0;for(j=2;j<=n;j+) bj=aj-u*bj-1-v*bj-2;r0=bn-1;r1=bn+u*bn-1;c0=b0;c1=b1-u*b0;for(j=2;j<=n-2;j+)cj=bj-u*cj-1-v*cj-2;s0=cn-3;s1=cn-2+u*cn-3;fru0=u*s0-s1;fru1=v*s0;frv0=-s0;frv1=-s1;deltau=-(frv1*r0-frv0*r1)/(fru0

19、*frv1-fru1*frv0);deltav=-(fru1*r0-fru0*r1)/(frv0*fru1-frv1*fru0);u=u+deltau;v+=deltav;delta=u*u-4*v;if(fabs(deltau)<eps)&&(fabs(deltav)<eps)print=1;printf("found roots:");r=-u/2;i=(int) sqrt(fabs(delta)/2;if(delta<0) printf("%.3f+%.3fit",r,fabs(double)i);printf(&

20、quot;%.3f-%.3fin",r,fabs(double)i);else printf("%.3ftt",r+i);printf("%.3fn",r-i);if(n=1) printf("Found root :%.3fn",a1/a0);if(k=2) return 0;if(k=1&&print) break; /* end while */ /* end for */return 0;运行:Please input max exponent:4Please input epslion:1e-5Pl

21、ease input the coefficient:1 1 5 4 4+1.000X4+1.000X3+5.000X2+4.000X1+4.000X0found roots:-0.500+0.000i-0.500-0.000ifound roots:0.000+2.000i0.000-2.000i习题3-15. 请编程用列全主元高斯约当消去法求矩阵A, B的逆矩阵。 #include <cstdio>#include <cmath>#define MAX 10int ROW;static float AMAX+12*MAX+1;void putout()int i,j

22、;printf("nThe reverse matrix is:n");for(i=1;i<=ROW;i+) for(j=ROW+1;j<=2*ROW;j+) printf("%.3f ",Aij);printf("n");int get()float maxMAX+1;int count_rMAX+1;float tmp2*MAX+1;float pass1,pass2;int i,j,k;for(i=1;i<=ROW;i+) maxi=Aii;count_ri=i;for(k=i;k<=ROW;k+) if

23、(maxi<Aki) maxi=Aki;count_ri=k; /* find the max one*/if(maxi=0) return -1;if(count_ri!=i) for(k=1;k<=2*ROW;k+)tmpk=Acount_rik;Acount_rik=Aik;Aik=tmpk; /*change the rows*/pass1=Aii;for(k=1;k<=2*ROW;k+) Aik/=pass1;for(j=1;j<=ROW;j+) if(j=i) continue;else pass2=Aji;for(k=1;k<=2*ROW;k+) A

24、jk=Ajk-pass2 *Aik; /* get simple*/putout();return 0;int main()int i,j;float tmp=-1;float p;printf("please input the number of ROW:");scanf("%d",&ROW);printf("n");if(ROW>=MAX) printf("the number of row is too big!n");return 0;printf("Please input th

25、e matrix A:n");for(i=1;i<=ROW;i+) for(j=1;j<=ROW;j+) scanf("%f,",&p);Aij=p;if(tmp<fabs(Aij) tmp=fabs(Aij);if(tmp=0) printf("No inverse!");return 0;for(i=1;i<=ROW;i+) Aii+ROW=1;return get();运行:please input the number of ROW:3Please input the matrix A:-3 8 52 -

26、7 41 9 -6The reverse matrix is:0.026 0.396 0.285 0.068 0.055 0.094 0.106 0.149 0.021please input the number of ROW:4Please input the matrix A:2 1 -3 -13 1 0 7-1 2 4 -21 0 -1 5The reverse matrix is:-0.047 0.588 -0.271 -0.9410.388 -0.353 0.482 0.765-0.224 0.294 -0.035 -0.471-0.035 -0.059 0.047 0.294习题

27、3-24. 编程用追赶法解下列三对角线方程组。#include <stdio.h>#include <math.h>#define ROW 5static float AROW+1=0,0,0,0,0,0;static float BROW+1=0,0,0,0,0,0;static float CROW+1=0,0,0,0,0,0;static float XROW+1=0,0,0,0,0,0;int main()int i;printf("Please input the A:n");for(i=2;i<=ROW;i+) scanf(&quo

28、t;%f,",&Ai);printf("Please input the B:n");for(i=1;i<=ROW;i+) scanf("%f,",&Bi);printf("Please input the C:n");for(i=1;i<=ROW-1;i+) scanf("%f,",&Ci);printf("Please input the F:n");for(i=1;i<=ROW;i+) scanf("%f,",&

29、;Xi);C1=C1/B1;for(i=2;i<=ROW-1;i+)Ci=Ci/(Bi-Ai*Ci-1);X1=X1/B1;for (i=2; i<=ROW; i+)Xi=(Xi-Ai*Xi-1)/(Bi-Ai*Ci-1);for(i=ROW-1;i>=1;i-)Xi=Xi-Ci*Xi+1;printf ("The value is : n");for(i=1;i<=ROW;i+) printf ("x%d=%.3fn",i,Xi);return 0;运行:Please input the A:-1 -1 -1 -1Please

30、input the B:4 4 4 4 4 Please input the C:-1 -1 -1 -1Please input the F:100 200 200 200 100The value is : x1=46.154x2=84.615x3=92.308x4=84.615x5=46.154习题4-31. 用高斯赛德尔方法编程解下列线性方程组,要求当时迭代终止。#include<cstdio>#include<cmath>float x6;float y6;float eps;float a66=4,-1,0,-1,0,0,-1,4,-1,0,-1,0,0,-1

31、,4,0,0,-1,-1,0,0,4,-1,0,0,-1,0,-1,4,-1,0,0,-1,0,-1,4;float b6=0,5,0,6,-2,6;int gs()int i,k,j;float s;for(i=0;i<6;i+) yi=xi=0;for(k=0;k<20;k+) for(i=0;i<6;i+) s=0;for(j=0;j<6;j+) if(j!=i) s=s+aij*yj;yi=(bi-s)/aii;for(i=0;(i<6)&&(abs(yi-xi)<eps);i+) ;if(i>=5) break;for(j=0

32、;j<6;j+) xj=yj;if(k>20) printf("No convergence./n");return -1;return 1;int main()int tag,i,m;printf("nthe max circles=");scanf("%d",&m);printf("eps=");scanf("%lf",&eps);tag=gs();if(tag>0) for(i=0;i<=5;i+) xi=yi;printf("x(%d)=

33、%.3en",i+1,xi);return 0;运行:nthe max circles=10eps=1e-5x(1)=1.000e+000x(2)=2.000e+000x(3)=1.000e+000x(4)=2.000e+000x(5)=1.000e+000x(6)=2.000e+000习题4-4#include <cmath>#include <cstdio>int agsdl(double a10, double b,int n, double x,double eps,double w,int m) int i,j,pcount; double p,t,

34、s,q; for(i=0;i<10;i+) xi=0; p=eps+1.0; for (pcount=0;p>=eps&&pcount<m;pcount+) p=0.0;for (i=0; i<n; i+) t=xi; s=0.0;for (j=0; j<n; j+)if (j!=i) s=s+aij*xj;xi=(bi-s)/aii;q=w*fabs(xi-t);if (q>p) p=q; return pcount;int main() int i,j,n,m,count=0;double eps,w3;static double a10

35、10;static double x10,b10;printf("nthe max circles=");scanf("%d",&m);printf("n=");scanf("%d",&n);printf("nmatrix A:n"); for(i=0;i<n;i+)for(j=0;j<n;j+) scanf("%lf",&aij); printf("nvector b:n");for(i=0;i<n;i+) sc

36、anf("%lf",&bi);printf("eps=");scanf("%lf",&eps);for(i=0;i<3;i+) printf("w(%d)=",i); scanf("%lf",&wi);for(j=0;j<3;j+) printf("w(%d)=%.3fn",j,wj);if(agsdl(a,b,n,x,eps,wj,m)<m) for(i=0;i<n;i+) printf("x(%d)=%.3e &q

37、uot;,i,xi);printf("n");elseprintf("failn");return 0;运行:the max circles=10n=3matrix A:4 -1 0-1 4 -10 -1 4vector b:1 4 -3eps=1e-5w(0)=1.03w(1)=1w(2)=1.1w(0)=1.030x(0)=5.000e-001 x(1)=1.000e+000 x(2)=-5.000e-001 w(1)=1.000x(0)=5.000e-001 x(1)=1.000e+000 x(2)=-5.000e-001 w(2)=1.100x(

38、0)=5.000e-001 x(1)=1.000e+000 x(2)=-5.000e-001 第5章习题5-12. 编程求下列矛盾方程组的解:#include<fstream.h>#include<stdio.h>#include<stdlib.h>double eps=1.0e-5;/error rangedouble *a; /save matrix Adouble *ab;/save augmented matrix ABdouble *ta;/save matrix ATAdouble *b; /save vectors Bdouble *tb; /

39、save ATBdouble *ca;/save matrix Adouble *cb; /save vectors Bdouble *x; /save result vectors Xint *fr;/row flagint *fc;/column flagint m,cm;/row num of matrix Aint n;/column num of matrix Aint r;void init();int G_J();int d_0(int);void change();void g_a_m();void show();void show1();void main()int t;in

40、it();t=1;while(1)r=G_J();if (m=r) break;if (d_0(r)=0) break;change();t=0;if (n=r)if (t)cout<<"该线性方程组是非矛盾线性方程组,有唯一解:"<<endl;show();elsecout<<"该线性方程组是矛盾线性方程组,有唯一最小二乘解:"<<endl;show();elseif (t)cout<<"该线性方程组是非矛盾线性方程组,有普遍解:"<<endl;show1();

41、elsecout<<"该线性方程组是矛盾线性方程组,有普遍最小二乘解:"<<endl;show1();void init() /open the data file-ab.txtifstream fin("ab.txt");if(!fin)cout<<"Can't open file ab.txt to read!"<<endl;exit(-1);/read the row and column information of matrix Afin>>m>>

42、;n;cm=m;/require room of matrix Aa=(double *)malloc(m+1)*sizeof(double *);int i,j;for(i=1;i<=m;i+)ai=(double *)malloc(n+1)*sizeof(double);/read matrix A's data from ab.txtfor(i=1;i<=m;i+)for(j=1;j<=n;j+)fin>>aij;/init matrix Bb=(double*)malloc(m+1)*sizeof(double);for(i=1;i<=m;i

43、+)fin>>bi;/close the data filefin.close();/backup matrix A and Bca=a; cb=b;return;/generate augmented matrixvoid g_a_m()int i,j;/require roomab=(double *)malloc(m+1)*sizeof(double*);for(i=1;i<=m;i+)abi=(double *)malloc(n+2)*sizeof(double);/load datafor(i=1;i<=m;i+)for(j=1;j<=n;j+)abij

44、=aij;abin+1=bi;/init the flag array fr(flag row) and fc(flag column)fr=(int*)malloc(m+1)*sizeof(int);fc=(int*)malloc(n+1)*sizeof(int);for(i=1;i<=m;i+) fri=i;for(j=1;j<=n;j+) fcj=j;return;/Gauss-Jordon elimination methodint G_J()int k,found,min;int i,j,pi,pj,ft;double t1,t2;/generate augmented

45、matrixg_a_m();k=1; min=(m<n)?m:n;/look for main elementwhile(k<=min)t1=found=pi=pj=0;for(i=k;i<=m;i+)for(j=k;j<=n;j+)/t2=abs(abfrifcj);t2=abfrifcj;t2=(t2>0)?t2:(-t2);if(t2>t1)&&(t2>eps)found=1; t1=t2;pi=i; pj=j;/end of j/end of iif(!found)return k-1;/found a main element

46、/fix the flag array-exchange the row and columnft=frpi; frpi=frk; frk=ft;ft=fcpj; fcpj=fck; fck=ft;/row standardizefor(j=1;j<=n+1;j+)abfrkj/=t1;/eliminationfor(i=1;i<=m;i+)if(i=k) continue;t2=abfrifck;for(j=1;j<=n+1;j+)abfrij-=abfrkj*t2;/look for next main elementk+;/end of whilereturn k-1;

47、/a=aTa, b=aTbvoid change()/require room for aTata=(double*)malloc(n+1)*sizeof(double*);int i;for(i=1;i<=n;i+)tai=(double*)malloc(n+1)*sizeof(double);/computing aTaint j,k;double sum;for(i=1;i<=n;i+)for(j=1;j<=n;j+)sum=0;for(k=1;k<=m;k+)sum+=aki*akj;taij=sum;/require roo for aTbtb=(double

48、*)malloc(n+1)*sizeof(double);/computing aTbfor(i=1;i<=n;i+)sum=0;for(k=1;k<=m;k+)sum+=aki*bk;tbi=sum;/free a,b/a=aTa, b=aTba=ta; b=tb;m=n;return;/if from d(r+1) to d(m) are all zero, return 0/otherwise,return 1int d_0(int r)int j;int t0;float t1;t0=0;for(j=r+1;j<=m;j+)t1=abfrjn+1;if(t1>e

49、ps|t1<-eps)t0=1;return t0;/show the augmented matrix ABvoid show()int i,j;/output and store the result X0.n-1x=(double*)malloc(n*sizeof(double);for(j=1;j<=n;j+)for(i=1;i<=m;i+)if(abij=1)xj-1=abin+1;cout<<"X"<<j-1<<"="<<xj-1<<endl;/output the

50、 residual vectorscout<<"该方程组的剩余向量是:"<<endl;cout<<"r=("double sum;for(i=1;i<=cm;i+)sum=0;for(j=1;j<=n;j+)sum+=caij*xj-1;cout<<cbi-sum;if(i!=cm) cout<<","cout<<")T"<<endl;return;void show1()int i,j,k,found;for(i=1;

51、i<=r;i+)cout<<"X"<<fci<<"="<<abfrin+1;for(j=1;j<=n-r;j+)cout<<"-X"<<j<<"*("<<abfrifcr+j<<")"cout<<endl;for(i=r+1;i<=n;i+)cout<<"X"<<fci<<"=X"<<fci<<""<<endl;cout<<"是否需要特解?(y/n)"char ch;cin>>ch;if(ch='n')|(ch='N') return;cout<<"请输入参数:"<<endl;/input variablesx=(double*)mal

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

最新文档

评论

0/150

提交评论