数值分析编程题c语言
- 格式:doc
- 大小:103.50 KB
- 文档页数:12
数值分析上机题Chapter 1(1)从大到小顺序#include<iostream>#include<math.h>#include<string>using namespace std;void main(){int N;cout<<"Please input N:"<<endl;cin>>N;cout<<"The number you input is:N="<<N<<endl;float sn=0;for(int i=2;i<=N;++i)sn+=1.0/(i*i-1);cout<<"The result you want is:"<<endl<<"Sn="<<sn<<endl; }(2)从小到大的顺序#include<iostream>#include<math.h>#include<string>using namespace std;void main(){int N;cout<<"Please input N:"<<endl;cin>>N;cout<<"The number you input is:N="<<N<<endl;float sn=0;for(int i=N;i>1;--i)sn+=1.0/(i*i-1);cout<<"The result you want is:"<<endl<<"Sn="<<sn<<endl; }(3)结果:N 100 10000 1000000(1) 0.740049 0.749852 -14.2546(2) 0.74005 0.7499 -14.2551 Chapter 2(1) 通用程序:#include <iostream>#include <string>#include <math.h>using namespace std;double h;//允许误差double x[100];double x0;//初值double xl;//所求结果double f(double x) //原函数{double a;//函数值a=f(x);return a;}double f1(double x)//导函数{double b=0.01,c[100];d;for(int i=0; ;++i){c[i]=(f(x+b)-f(x))/b;//导函数if(i>=1){if(fabs(c[i]-c[i-1])<=h){d=c[i];break;}}b/=10;}return d;}void main(){for(int k=0; ;++i){x[0]=x0;//初值x[k+1]=x[k]-f(x[k])/f1(x[k]);if(fabs(x[k+1]-x[k])<=h)break;}cout<<"xl="<<x[k+1]<<endl;//输出结果}(2)#include <iostream>#include <string>#include <math.h>using namespace std;int k;double h=0.001;//允许误差double x[100];double x0;//初值double xl;//所求结果double f(double x) //原函数{double a;//函数值a=x*x*x/3-x;return a;}double f1(double x)//导函数{double b=0.01,c[100],d;for(int i=0; ;++i){c[i]=(f(x+b)-f(x))/b;//导函数if(i>=1){if(fabs(c[i]-c[i-1])<=h){d=c[i];break;}}b/=10;}return d;}int main(){for(int i=1;i<=1000;++i){x0=i/1000.0;for(k=0; ;++k){x[0]=x0;x[k+1]=x[k]-f(x[k])/f1(x[k]);if(fabs(x[k+1]-x[k])<=h)break;}if(fabs(x[k+1])>=1)break;}cout<<"when deta="<<x0<<": "<<"xl*="<<x[k+1]<<endl;//输出结果cout<<"This is the biggest deta,or the result will not be zero."<<endl;return 0;}Chapter 3(1).通用程序:#include <iostream>#include <string>#include <math.h>#define n 3//定义阶次using namespace std;void main(){int i,j;double a[n][n+1],temp,x[n];cout<<"请输入需求解矩阵:"<<endl;//输入求解矩阵for(i=0;i<n;++i)for(j=0;j<=n;++j)cin>>a[i][j];for(j=0;j<n;++j)//列主元{ int k=j;double h=a[j][j];//列首元int b=j;for(i=j+1;i<n;++i)//列最大元素并标记{if(fabs(a[i][k])>fabs(h)){h=a[i][k];b=i;}}for(int m=j;m<=n;++m)//交换主元{temp=a[k][m];a[k][m]=a[b][m];a[b][m]=temp;}for(int c=j+1;c<n;++c)//化为上三角矩阵for(int z=j;z<=n;++z){double yz=a[c][j]/h;a[c][z]=a[c][z]-a[j][z]*yz;}}//求解过程x[n-1]=a[n-1][n]/a[n-1][n-1];cout<<"x"<<"["<<n-1<<"]="<<x[n-1]<<endl;for(int e=n-2;e>=0;--e){ double l=0;for(int w=e+1;w<n;++w){l=l+a[e][w]*x[w];}x[e]=(a[e][n]-l)/a[e][e];cout<<"x"<<"["<<e<<"]="<<x[e]<<endl;}}x0=-0.23645,x1=0.47743,x2=-0.76698,x3=-0.077649,x4=-0.30792,x5=0.14634,x6=-0.17 073,x7=0.28480,x8=0.34483;Chapter 437.(1)通用程序程序在使用前要先设定n的值#include<iostream>#include<math.h>#include<string>using namespace std;#define n 10double x[n+1],y[n+1],h[n],u[n],u_u[n],d[n+1];double a[n+1][n+1],m[n+1];double F0,Fn,arr;double b[n+1][n+1];double f(double &a){double function;for(int i=0;i<n+1;++i)if(a==x[i])function=y[i];return function;}double f1(double &a,double &b){if(a!=b)return (f(b)-f(a))/(b-a);else{if(a==x[0])return F0;elsereturn Fn;}}double f2(double &a,double &b,double &c){return (f1(b,c)-f1(a,b))/(c-a);}double s(double x0){double sx;if(x0<=x[0]||x0>=x[n]){cout<<"不?在¨²定¡§义°?域®¨°内¨²"<<endl;return 0;}else{int j;for(int i=0;i<n+1;++i)if(x0>=x[i]&&x0<=x[i+1]) j=i;sx=y[j]+(f1(x[j],x[j+1])-(m[j]/3.0+m[j+1]/6.0)*h[j])*(x0-x[j])+(m[j]/2.0)*(x0-x[j])*(x0-x [j])+((m[j+1]-m[j])/6*h[j])*(x0-x[j])*(x0-x[j])*(x0-x[j]);return sx;}}int main(){cout<<"please input x:";for(int i1=0;i1<n+1;++i1)cin>>x[i1];cout<<endl<<endl;cout<<"please input y:";for(int i2=0;i2<n+1;++i2)cin>>y[i2];cout<<endl<<endl;cout<<"please input F0,Fn:";cin>>F0>>Fn;cout<<endl<<endl;for(int i3=0;i3<n;++i3)h[i3]=x[i3+1]-x[i3];for(int i4=1;i4<n;++i4){u[i4]=h[i4-1]/(h[i4-1]+h[i4]);u_u[i4]=1-u[i4];}for(int i=0;i<n+1;++i)for(int j=0;j<n+1;++j){if(i==j)a[i][j]=2;else if(j==i-1){if(i!=n)a[i][j]=u[i];elsea[i][j]=1;}else if(j==i+1){if(i!=0){a[i][j]=u_u[i];}elsea[i][j]=1;}elsea[i][j]=0;cout<<a[i][j]<<" ";if(j==n)cout<<endl<<endl;}cout<<"d[i]的Ì?值¦Ì:êo"<<endl;for(int i5=0;i5<n+1;++i5){if(i5==0)d[i5]=6*f2(x[0],x[0],x[1]);else if(i5==n)d[i5]=6*f2(x[n-1],x[n],x[n]);elsed[i5]=6*f2(x[i5-1],x[i5],x[i5+1]); cout<<d[i5]<<" ";}cout<<endl<<endl;cout<<"系¦Ì数ºy矩?阵¨®:êo"<<endl<<endl; for(int i=0;i<n+1;++i)for(int j=0;j<=n+1;++j){if(j<n+1)b[i][j]=a[i][j];elseb[i][j]=d[i];cout<<b[i][j]<<" ";if(j==n+1)cout<<endl;}cout<<endl<<endl;for (int i=1;i<=n;i++){ arr=b[i][i-1];for (int j=0;j<=n+1;j++)b[i][j]=b[i][j]-a[i-1][j]*arr/b[i-1][i-1]; }m[n]=b[n][n+1]/b[n][n];for(int i=n-1;i>=0;i--) m[i]=(b[i][n+1]-b[i][i+1]*m[i+1])/b[i][i];for(int i=0;i<n+1;++i)cout<<m[i]<<" ";cout<<endl<<endl;for(int j=0;j<n;++j){cout<<"当Ì¡Àx在¨²区?间?["<<x[j]<<","<<x[j+1]<<"]时º¡À:";cout<<"S(x)="<<y[j]<<"+"<<(f1(x[j],x[j+1])-(m[j]/3.0+m[j+1]/6.0)*h[j]) <<"(x-"<<x[j]<<")+"<<m[j]/2.0<<"(x-"<<x[j]<<")^2+"<<(m[j+1]-m[j])/6*h[j]<<"(x-"<<x[j]<<")^3"<<endl<<endl;}for(int i=0;i<10;++i)cout<<"s("<<i<<"+0.5)="<<s(i+0.5)<<" "<<endl;return 0;}(2)此时n为10,将程序中的n设为10,输入数据运行,得到结果:Chapter 523(1)通用程序程序在运行之前要先设定积分区间,即a,b,c,d的值,本题为了编程方便,已经把数据写入程序中。
第一章作业#include<stdio.h>const int N=100;int main(){int n;float S1,S2;double S3;S1=0;S2=0;for(n=2;n<N+1;n++){S1=S1+1.0/(n*n-1);}printf("the sum S1 is %d\n",S1);for(n=N;n>1;n--){S2=S2+1.0/(n*n-1);}printf("the sum S2 is %d\n",S2);S3=0.5*(1.5-1.0/N-1.0/(N+1));printf("S3 is %d\n",S3);}修改N的值即可第二章作业(1)float df(float x){float df;df=x*x-1;return df;}main(){float x0,x1,a;int k=0;printf("请输入初值x0=");scanf_s("%f",&x0);do{a=-f(x0)/df(x0);x1=x0+a;k++;x0=x1;}while(fabs(a)>eps);printf("迭代次数k=%d,根值x0=%f",k,x0);}(2)#include<stdio.h>#include<math.h>#define eps 0.01float f(float x){float f;f=x*x*x/3-x;return f;}float df(float x){float df;df=x*x-1;return df;}float ddf(float x){float ddf;ddf=2*x;return ddf;}main(){float x0,x1,a,dt=0;int k=0;do{dt=dt+eps;k++;}while((f(-dt)*f(dt)<0)&&(df(dt)!=0)&&(2*dt>=-f(-dt)/df(-dt))&&(2*dt>=f(dt)/df(dt)));printf("请输入x0:");scanf_s("%f",&x0);if(x0<-1)do{a=-f(x0)/df(x0);x1=x0+a;k++;x0=x1;}while(fabs(a)>eps);else if(x0>-1&&x0<-dt)do{a=-f(x0)/df(x0);x1=x0+a;k++;x0=x1;}while(fabs(a)>eps);else if(x0>-dt&&x0<dt)do{a=-f(x0)/df(x0);x1=x0+a;k++;x0=x1;}while(fabs(a)>eps);else if(x0>dt&&x0<1)do{a=-f(x0)/df(x0);x1=x0+a;k++;x0=x1;}while(fabs(a)>eps);else if (x0>1)do{a=-f(x0)/df(x0);x1=x0+a;k++;x0=x1;}while(fabs(a)>eps);printf("收敛域=%f,迭代次数=%d\n",dt,k);printf("根=%f\n",x0);}第三章作业#include<stdio.h>#include<math.h>main(){int n,i,j,k,t;float a[10][11],s[10],x[10],s2,max,sum,b;printf("请输入n= ");scanf("%d",&n);printf("请输入矩阵A= \n");for(i=0;i<=n;i++)for(j=0;j<=(n+1);j++)scanf("%f",&a[i][j]);//变换成高斯消去法矩阵//for(k=0;k<n;k++){//找出第k列中绝对值最大的一项//j=0;for(i=k;i<=n;i++){s[j]=a[i][k];j++;}max=fabs(s[0]);t=k;for(i=1;i<j;i++){if (fabs(s[i])>max){max=fabs(s[i]);t=k+i;}}// printf("t=%d\n",t);//将在列中有最大的元素的一行和第k行交换//for(j=0;j<=n+1;j++){b=a[k][j];a[k][j]=a[t][j];a[t][j]=b;}/*printf("***\n");for(i=0; i<=n; i++){for(j=0; j<=(n+1); j++)printf("%f ",a[i][j]);printf("\n");}*///将第k行k列的元素变为零//for(i=(k+1);i<=n;i++){s2=a[i][k]/a[k][k];for(j=k;j<=(n+1);j++)a[i][j]=a[i][j]-s2*a[k][j];}/*printf("###\n");for(i=0; i<=n; i++){for(j=0; j<=(n+1); j++)printf("%f ",a[i][j]);printf("\n");}*/printf("\n");}printf("\n");//解方程--//for(i=n;i>=0;i--){sum=0;x[n+1]=a[n][n+1]/a[n][n];for(j=i;j<=n;j++)sum+=a[i-1][j]*x[j+1];x[i]=(a[i-1][n+1]-sum)/a[i-1][i-1];}printf("方程的解为:");for(i=0;i<=n;i++)printf("%f,",x[i+1]);}。
上机实习题一一、题目:已知A 与b12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.0317432.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 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,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.6784392} 1.用Household 变换,把A 化为三对角阵B (并打印B )。
C语言编程题带答案题目 1:求两个整数的最大值```cinclude <stdioh>int max(int num1, int num2) {if (num1 > num2) {return num1;} else {return num2;}}int main(){int num1 = 10, num2 = 20;int maxValue = max(num1, num2);printf("最大值为: %d\n", maxValue);return 0;}```分析:在这个程序中,我们定义了一个名为`max` 的函数,它接受两个整数参数`num1` 和`num2` 。
通过使用条件判断语句`if` 来比较这两个数的大小,如果`num1` 大于`num2` ,则返回`num1` ,否则返回`num2` 。
在`main` 函数中,我们给定了两个整数`num1` 和`num2` 的值,并调用`max` 函数来获取它们中的最大值,最后使用`printf` 函数将最大值输出到控制台。
题目 2:计算一个整数数组的平均值```cinclude <stdioh>float average(int arr, int size) {int sum = 0;for (int i = 0; i < size; i++){sum += arri;}return (float)sum / size;}int main(){int arr ={10, 20, 30, 40, 50};int size = sizeof(arr) / sizeof(arr0);float avg = average(arr, size);printf("平均值为: %2f\n", avg);return 0;}```分析:在这个程序中,首先在`average` 函数里,我们初始化一个变量`sum` 为 0 ,用于存储数组元素的总和。
《数值分析》计算实习作业《一》北航第一题 设有501501⨯的矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=501500499321a bc b a b cc b a b ccb a bc c b a b c b a A其中.064.0,16.0);501,2,1(64.0)2.0sin()024.064.1(1.0-===--=c b i e i i a i i 矩阵的特征值)501,,2,1( =i i λ满足||min ||,501150121i i s λλλλλ≤≤=<<<试求1. 5011,λλ和s λ的值2. 的与数4015011λλκλμ-+=k 最接近的特征值)39,,2,1( =K κλi3. 的(谱范数)条件数2)A (cond 和行列式A det 要求1. 算法的设计方案(A 的所有零元素都不能存储)2. 全部源程序(详细注释)。
变量为double ,精度-1210=ε,输出为e 型12位有效数字3. 特征值s 5011,,λλλ和)39,,2,1( =K κλi 以及A cond det ,)A (2的值4. 讨论迭代初始向量的选取对计算结果的影响,并说明原因解答:1. 算法设计对于s λ满足||min ||5011i i s λλ≤≤=,所以s λ是按模最小的特征值,直接运用反幂法可求得。
对于5011,λλ,一个是最大的特征值,一个是最小的特征值,不能确定两者的绝对值是否相等,因此必须首先假设||||5011λλ≠,然后运用幂法,看能否求得一个特征值,如果可以求得一个,证明A 是收敛的,求得的结果是正确的,然后对A 进行带原点平移的幂法,偏移量是前面求得的特征值,可以求得另一个特征值,最后比较这两个特征值,较大的特征值是501λ,较小的特征值就是1λ。
如果在假设的前提下,无法运用幂法求得按模最大的特征值,即此时A 不收敛,则需要将A 进行带原点平移的幂法,平移量可以选取1,再重复上述步骤即可求得两个特征值。
拉格朗日插值#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h>void Lagra nge(float s){double x[5]={0.2,0.4,0.6,0.8,1.0},y[5]={0.98,0.92,0.81,0.64,0.38},f,L=0; int i,j;for (i=0;i<5;i++){ f=1; for (j=0;j<5;j++) if(j!=i) f=(s-x[j])/(x[i]-x[j])*f; L+=f*y[i]; } printf("输出:%f\n",L);}void mai n(){float x;printf("输入插值点:");scan f("%f", &x);Lagra nge(x);}妙人斑值原:0-5爺出:0-871406any key to continue牛顿插值#in clude<stdlib.h>#in clude<stdio.h>#in clude<math.h> int ND(float s){double x[ 5]={0.2,0.4,0.6,0.8,1.0},y[5]={0.98,0.92,0.81,0.64,0.38},p=0,g,f;int i,j,k;for (i=0;i<5;i++){for (j=4;j>i;j--){ f=x[j]-x[j-i-1];y[j]=(y[j]-y[j-1])/f;}g=y[i+1];for (k=0;k<=i;k++)g=g*(s-x[k]);p=p+g;}printf("输出插值点函数值:%f\n",p+y[0]);return 1;void mai n(){float x;printf("输入插值点:"); scan f("%f", &x);ND(x);}三、埃尔米特插值#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h>void Hermite(float s){double x[3]={0.25,1,2.25},y[3]={0.125,1,3.375},z[3]={0.75,1.5,2.25};double H=0.0,a,b,f,g;int i,j;for (i=0;i<3;i++){f=1.0;g=0.0;for (j=0;j<3&&j!=i;j++) {f=f*(s-x[j])/(x[i]-x[j]);g=g+1/(x[i]-x[j]);} a=(1-2*(s-x[i])*g)*f*f;b=(s-x[i])*f*f;H=H+y[i]*a+z[i]*b;}prin tf("%f\n",H);}void mai n(){float x;printf("输入插值点:");scan f("%f", &x);Hermite(x);}四、三次样条插值#i nclude <math.h>#i nclude <stdio.h>#in clude <stdlib.h>void main(){int N=7,R=2,i,k;double p1,p2,p3,p4;double x[8]={0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9};double y[8]={0.4794,0.6442,0.7833,0.8912,0.9636,0.9975,0.9917,0.9463};double P0=-0.4794,Pn=-0.9463,u[3]={0.6,0.8,1.2},s[3];double h[7],a[8],c[7], g[8],af[8],ba[7],m[8];for(k=0;k <N;k++)h[k]=x[k+1]-x[k];for(k=1;k <N;k++) for(k=1;k <N;k++) for(k=1;k <N;k++) c[0]=a[N]=1; a[k]=h[k]/(h[k]+h[k-1]);c[k]=1-a[k];g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]);g[0]=3*(y[1]-y[0])/h[0]-P0*h[0]/2;g[N]=3*(y[N]-y[N-1])/h[N-1]+Pn*h[N-1]/2;ba[0]=c[0]/2; g[0]=g[0]/2; for(i=1;i <N;i++){ af[i]=2-a[i]*ba[i-1]; g[i]=(g[i]-a[i]*g[i-1])/af[i]; ba[i]=c[i]/af[i]; } af[N]=2-a[N]*ba[N-1]; g[N]=(g[N]-a[N]*g[N-1])/af[N]; m[N]=g[N];for(i=N-1;i>=0;i--) m[i]=g[i]-ba[i]*m[i+1];for(i=0;i <=R;i++){ k=0;while(u[i]> x[k+1])k++;p1=(h[k]+2*(u[i]-x[k])*pow((u[i]-x[k+1]),2)*y[k])/pow(h[k],3); p2=(h[k]-2*(u[i]-x[k+1])*pow((u[i]-x[k]),2)*y[k+1])/pow(h[k],3); p3=(u[i]-x[k])*pow((u[i]-x[k+1]),2)*m[k]/pow(h[k],2); p4=(u[i]-x[k+1])*pow((u[i]-x[k]),2)*m[k+1]/pow(h[k],2); s[i]=p1+p2+p3+p4;}printf( "\nx= ");for(i=0;i <=N;i++)printf( "%8.1f ",x[i]);printf( "\ny= ");for(i=0;i <=N;i++)printf( "%8.4f ",y[i]);printf( "\n\nu= ");for(i=0;i <=R;i++)printf( "%9.2f ",u[i]);printf( "\n 插值点:s= ");for(i=0;i <=R;i++)prin tf( "%9.5f ”,s[i]);prin tf("\n");五、复合梯形公式#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h> double FTX(i nt n,float a,float b){double f=0,t,h,*x,*y;int i;x=(double*)malloc(( n+1)*sizeof(double)); y=(double*)malloc(( n+1)*sizeof(double)); h=(b-a)/n;for(i=0;i< n+1;i++){x[i]=a+i*h;y[i]=s in (x[i]);}for(i=1;i< n; i++) f=f+2*y[i];t=h/2*(y[0]+f+y[ n]);printf("输出函数值:%f\n",t);return 1;}void mai n(){float a,b;int n;printf("输入区间上,下限:"); scan f("%f %f", &a, &b);printf("输入等分区间数:");scan f("%d",&n);FTX( n,a,b);}六、复合辛普森求积公式#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h> double FSP(int n,float a,float b){double f1=0,f2=0,h,*x1,*y1,*x2,*y2; int i;x1= (double*)malloc(( n+1)*sizeof(double));y1= (double*)malloc(( n+1)*sizeof(double));x2=(double*)malloc( n*sizeof(double));y2=(double*)malloc( n*sizeof(double)); h=(b-a)/n;for(i=0;i< n+1;i++) {x1[i]=a+h*i;y1[i]=si n(x1[i])*x1[i];} for(i=0;i< n;i++){x2[i]=x1[i]+h/2;y2[i]=si n( x2[i])*x2[i];}for(i=1;i< n;i++) f仁f1+2*y1[i];for(i=0;i< n;i++) f2=f2+4*y2[i];printf(” 输出函数值:%f\n",h/6*(y1[0]+f1+f2+y1[n]));return 1;}void mai n(){float a,b;int n;printf("输入区间上,下限:");scan f("%f %f", &a, &b);printf("输入等分区间数:");scan f("%d",&n);FSP( n,a,b);}七、直接三角分解法#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h>void mai n(){doubleA[3][3]={0.25,0.2,0.166667,0.3333,0.25,0.2,0.5,1,2},x[3],y[3],b[3]={9,8,8},L[3][3],U[3][3],f1=0, f2=0; int i,j,k;for(i=0;i<3;i++)for(j=0;j<3;j++){U[i][j]=O;L[i][j]=O;}for(i=0;i<3;i++){U[0][i]=A[0][i];L[i][0]=A[i][0]/U[0][0];L[i][i]=1;}for(i=1;i<3;i++)for(j=i;j<3;j++){for(k=0;k<=i-1;k++){f 1= f1+L[i][k]*U[k][j];f2+=L[j][k]*U[k][i];}U[i][j]=A[i][j]-f1;L[j][i]=(A[j][i]-f2)/U[i][i];f1=0;f2=0;}y[O]=b[O];for(i=1;i<3;i++){for(j=0;j<=i-1;j++) f1+=L[i][j]*y[j];y[i]=b[i]-f1;f1=0;}x[2]=y[2]/U[2][2];for(i=1;i>=0;i__){for(j=i+1;j<3;j++) f2+=U[i][j]*x[j];x[i]=(y[i]-f2)/U[i][i];f2=0;} printf(” 输出L 矩阵:\n”);for(i=0;i<3;i++){for(j=0;j<3;j++)prin tf("%f ",L[i][j]);prin tf("\n");}printf(” 输出U 矩阵:\n");for(i=0;i<3;i++){for(j=0;j<3;j++)prin tf("%f ",U[i][j]);prin tf("\n");}printf(”输出求解结果:\n”);for(i=0;i<3;i++)prin tf("%f ",x[i]);prin tf("\n");八、改进的平方法#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h>void mai n(){double A[3][3]={2,-1,1,-1,-2,3,1,3,1},x[3],y[3],b[3]={4,5,6},d[3 ],L [3][3],U[3][3],f1=0,f2=0;int i,j,k, n=3;for(i=0;i< n;i++)for(j=0;j< n;j++) {U[i][j]=0;L[i][j]=0;}d[0]=A[0][0];L[0][0]=1;for(i=1;i< n; i++){L[i][i]=1;for(j=0;j<=i-1;j++) {for(k=0;k<=j-1;k++) f1+=U[i][k]*L[j][k];U[i][j]=A[i][j]-f1;L[i][j]=U[i][j]/d[j];f2+=U[i][j]*L[i][j];f1=0;} d[i]=A[i][j]-f2;f2=0;} y[0]=b[0];for(i=1;i< n; i++) {for(j=0;j<=i-1;j++) f1+=L[i][j]*y[j];y[i]=b[i]-f1;f1=0;}x[ n-1]=y[ n-1]/d[ n-1];for(i=n-2;i>=0;i--) {for(j=i+1;j< n;j++) f2+=L[j][i]*x[j];x[i]=y[i]/d[i]-f2;f2=0;}printf(” 输出L 矩阵:\n");for(i=0;i< n;i++){for(j=0;j< n;j++) prin tf("%f ",L[i][j]); prin tf("\n");}printf("输出U 矩阵:\n");for(i=0;i< n;i++){for(j=0;j< n;j++) prin tf("%f ",U[i][j]); prin tf("\n");}printf("输出求解结果:\n");for(i=0;i< n;i++) prin tf("%f ",x[i]);prin tf("\n");}九、追赶法double A[5][5]={2,-1,0,0,0,-1,2,-1,0,0,0,-1,2,-1,0,0,0,-1,2,-1,0,0,0,-1,2},f[5]={1,0,0,0,0};double x[5],y[5],U[5][5];int i,j, n=5;for(i=0;i< n;i++)for(j=0;j<n;j++) U[i][j]=0;U[0][0]=U[ n-1][ n-1]=1;U[0][1]=A[0][1]/A[0][0];for(i=1;i< n-1;i++){ U[i][i]=1;U[i][i+1]=A[i][i+1]/(A[i][i]-A[i][i-1]*U[i-1][i]);}y[0]=f[0]/A[0][0];for(i=1;i< n; i++)y[i]=(f[i]-A[i][i-1]*y[i-1])/(A[i][i]-A[i][i-1]*U[i-1][i]);x[ n_1]=y[ n-1];for(i=n-2;i>=0;i--) x[i]=y[i]-U[i][i+1]*x[i+1];printf("输出U 矩阵:\n");for(i=0;i< n;i++){ for(j=0;j< n;j++) prin tf("%f ",U[i][j]); prin tf("\n");}printf(”输出求解结果:\n");for(i=0;i< n;i++) prin tf("%f ",x[i]);prin tf("\n");}十、雅可比迭代法#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h>void mai n(){double A[3][3]={5,2,1,-1,4,2,2,-3,10},x[50][3],b[3]={-12,20,3},f=0,L=1;int n=3,i,j,k=1;for (i=0;i<3;i++) x[0][i]=0;/* 初始值*/while(L<0.003){for(i=0;i <n ;i++){ for(j=0;j< n;j++)if(j!=i) f=f+A[i][j]*x[k-1][j];x[k][i]=(b[i]-f)/A[i][i];}L=x[k][0]-x[k-1][0];for(i=1;i <n ;i++)if((x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])>1||(x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])<-1) L=x[k][i]-x[k-1][i];k++;}printf(”输出每次迭代求得的X值:\n");for(i=1;i<k;i++){printf("第%d 次迭代:”,i);for(j=0;j< n;j++)prin tf("%f ",x[i][j]);prin tf("\n ”);}printf("\n 输出迭代次数:%d\n",k-1);咼斯一塞德尔迭代法double A[3][3]={5,2,1,-1,4,2,2,-3,10},x[50][3],b[3]={-12,20,3},f1=0,f2=0, L=1;int i,j,k=1, n=3;for (i=0;i<3;i++) x[0][i]=0;while(L>0.00001||L<-0.00001){for(i=0;i< n;i++){ for(j=0;j< n;j++){if(j<i) f1+=A[i][j]*x[k][j];else if(j>i) f2+=A[i][j]*x[k-1][j];}x[k][i]=(b[i]-f1-f2)/A[i][i];}L=x[k][0]-x[k-1][0];for(i=1;i <n ;i++)if((x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])>1||(x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])<-1) L=x[k][i]-x[k-1][i];k++;}printf(”输出迭代X矩阵:\n"); for(i=1;i<k;i++){for(j=0;j< n;j++) prin tf("%f ",x[i][j]);prin tf("\n");}printf("\n 输出迭代次数:%d\n",k-1);double A[3][3]={5,2,1,-1,4,2,2,-3,10},x[30][3],b[3]={-12,20,3},f1=0,f2=0, L=1,w=0.9; inti,j,k=1, n=3;for(i=0;i<3;i++) x[0][i]=0;while(L>0.0001||L<-0.0001){for(i=0;i <n ;i++){ for(j=0;j< n;j++){if(j<i) f1+=A[i][j]*x[k][j];else if(j>i) f2+=A[i][j]*x[k-1][j];}x[k][i]=w*(b[i]-f1-f2)/A[i][i];}L=x[k][0]-x[k-1][0];for(i=1;i< n; i++) if((x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])>1||(x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])<-1)L=x[k][i]-x[k-1][i];k++;}printf(”输出迭代X矩阵:\n");for(i=1;i<k;i++){for(j=0;j< n;j++) prin tf("%f ",x[i][j]); prin tf("\n");}printf("\n 输出迭代次数:%d\n",k-1);}数值分析算法程序班级:计算08Q2 班姓名:甄彦福。
#include <stdio.h>#include <math.h>double f(double x){double ans;ans=exp(x);return ans;}void main(){double a=1,b=3,error=0.0001,t[20][20],h,c;int i,j,k,m,n;h=b-a;t[0][0]=h*(f(a)+f(b))/2;k=1;while(1){t[0][k]=0;m=1;for(j=0;j<k-1;j++)m=m*2;for(i=1;i<=m;i++)t[0][k]=t[0][k]+h*f(a+(i-0.5)*h);t[0][k]=(t[0][k]+t[0][k-1])/2;for(j=1;j<=k;j++){ c=1;for(n=0;n<j;n++)c=c*4;t[j][k-j]=(c*t[j-1][k-j+1]-t[j-1][k-j])/(c-1);}if(fabs(t[k][0]-t[k-1][0])<error){ printf("\n积分结果I ≈%lf\n",t[k][0]);break;}else{ h=h/2;k++;}}}#include <stdio.h>#include <math.h>double f(double t){double ans;ans=pow(cos(t),1.0/3);return ans;}void main(){double x=0,eslong=0.000001,x0;int N=20,i;printf("\n近似初值x0 = %lf\n",x);for(i=0;i<N;i++){x0=x;x=f(x);printf(" x%d = %lf\n",i+1,x);if(fabs(x-x0)<eslong)break;}if(fabs(x-x0)<eslong)printf("得到近似结果为x ≈%lf\n\n",x,i);elseprintf("迭代失败\n");}#include <stdio.h>#include <math.h>double a=0,b=1,x,y=0,h=0.1,k1,k2,k3,k4;int i,N;double f(double t,double s){double ans;ans=1+t*t;return ans;}void main(){N=(b-a)/h;x=a;printf("\n 初值为(x0,y0) = ( %.8f , %.8f )\n",x,y);for(i=0;i<N;i++){k1=f(x,y);k2=f(x+h/2,y+h*k1/2);k3=f(x+h/2,y+h*k2/2);k4=f(x+h,y+h*k3);y=y+h*(k1+2*(k2+k3)+k4)/6;x=x+h;printf(" 第%d次输出结果为(x%d,y%d) = ( %.8f , %.8f )\n",i+1,i+1,i+1,x,y);}}#include <stdio.h>void main(){double datax[4]={1.2,2.9,4.6,5.8},datay[10]={14.84,33.71,58.36,79.24},l[3],x=1.5,y;int i,j;y=0;for(i=0;i<=3;i++){l[i]=1;for(j=0;j<i;j++)l[i]=(x-datax[j])/(datax[i]-datax[j])*l[i];for(j=i+1;j<=3;j++)l[i]=(x-datax[j])/(datax[i]-datax[j])*l[i];y=y+datay[i]*l[i];}printf("\n f(x)在x = %f 处的近似值为: y = %f\n",x,y);}#include <stdio.h>void main(){double datay[9]={11.7,14.87,21.44,31.39,44.73,61.46,81.57,105.11,131.91};int m=2,i,j,k;double p,data[9][4],a[3][4],datax[9]={1.2,2.3,3.4,4.5,5.6,6.7,7.8,8.9,10.0};for(i=0;i<9;i++)for(j=1;j<2*m+1;j++){data[i][j]=1;for(k=0;k<j;k++)data[i][j]=data[i][j]*datax[i];}for(i=0;i<m+1;i++){for(j=0;j<m+1;j++){a[i][j]=0;for(k=0;k<9;k++)a[i][j]=a[i][j]+data[k][i+j];}}a[0][0]=9;a[0][m+1]=0;for(i=0;i<9;i++)a[0][m+1]=a[0][m+1]+datay[i];for(i=1;i<m+1;i++){ a[i][m+1]=0;for(j=0;j<9;j++){p=datay[j];for(k=0;k<i;k++)p=p*datax[j];a[i][m+1]=a[i][m+1]+p;}} //生成m+1行,m+2列增广矩阵// for(i=0;i<m+1;i++) //显示方程组//for(j=0;j<m+2;j++){if(j!=m+1){ printf("(%f)a%d ",a[i][j],j);if(j!=m)printf("+ ");}elseprintf("= %f \n",a[i][j]);}for(i=0;i<m;i++) //高斯消去法//{if(a[i][i]!=0){ for(j=i+1;j<m+1;j++){ a[j][i]=a[j][i]/a[i][i];for(k=i+1;k<m+2;k++)a[j][k]=a[j][k]-a[i][k]*a[j][i];}}elsebreak;}if(a[m][m]!=0&&i==m){ a[m][m+1]=a[m][m+1]/a[m][m];for(i=2;i<=m+1;i++){ for(j=1;j<i;j++)a[m+1-i][m+1]=a[m+1-i][m+1]-a[m+1-i][m+1-j]*a[m+1-j][m+1];a[m+1-i][m+1]=a[m+1-i][m+1]/a[m+1-i][m+1-i];}printf("方程组的解为:\n");for(j=0;j<m+1;j++)printf("a%d = %f\n",j,a[j][m+1]);printf("拟合多项式为:\n");printf("P%d(x) = (%f) + (%f)x + (%f)x^2\n",m,a[0][m+1],a[1][m+1],a[2][m+1]);}elseprintf("数据有误!\n");}}列主元素法#include <stdio.h>#include <math.h>void main(){double a[3][4]={1,-2,-1,3,-2,10,-3,15,-1,-2,5,10},mov,comp;int i,j,k,nrow;for(i=0;i<2;i++){comp=fabs(a[i][i]);for(k=i;k<3;k++) //比较绝对值大小并进行主元列交换// if(fabs(a[k][i])>=comp){nrow=k;comp=fabs(a[k][i]);}for(j=0;j<=3;j++){mov=a[i][j];a[i][j]=a[nrow][j];a[nrow][j]=mov;}printf("方程第%d行互换位置后如下\n",i+1);for(j=0;j<3;j++)printf("(%f)x1 + (%f)x2 + (%f)x3 = %f\n",a[j][0],a[j][1],a[j][2],a[j][3]);if(a[i][i]!=0){for(j=i+1;j<3;j++){a[j][i]=a[j][i]/a[i][i];for(k=i+1;k<=3;k++)a[j][k]=a[j][k]-a[i][k]*a[j][i];a[j][i]=0;}printf("方程经%d次消元如下\n",i+1);for(j=0;j<3;j++)printf("(%f)x1 + (%f)x2 + (%f)x3 = %f\n",a[j][0],a[j][1],a[j][2],a[j][3]);}elsebreak;}if(a[2][2]!=0&&i==2){printf("方程化简得\n");for(i=0;i<3;i++)printf("(%f)x1 + (%f)x2 + (%f)x3 = %f\n",a[i][0],a[i][1],a[i][2],a[i][3]);a[2][3]=a[2][3]/a[2][2];for(i=2;i<=3;i++){for(j=1;j<i;j++)a[3-i][3]=a[3-i][3]-a[3-i][3-j]*a[3-j][3];a[3-i][3]=a[3-i][3]/a[3-i][3-i];}printf("方程组的解为:\n");for(j=0;j<3;j++)printf("x%d = %f\n",j+1,a[j][3]);}elseprintf("数据有误!\n");}Jacobi迭代法#include <stdio.h>#include <math.h>void main(){double a[3][7]={{1,-2,-1,3},{-2,10,-3,15},{-1,-2,5,10}},error=0.000001,norm;int N=423,i,j,k;a[0][4]=0,a[1][4]=0,a[2][4]=0;for(i=0;i<3;i++) //把a矩阵转化为b矩阵//{a[i][6]=a[i][i];for(j=0;j<3;j++){a[i][j]=-a[i][j]/a[i][6];}a[i][3]=a[i][3]/a[i][6];a[i][i]=0;}printf("化为b矩阵如下\n");for(i=0;i<3;i++){printf("%f %f %f %f\n",a[i][0],a[i][1],a[i][2],a[i][3]);}for(i=1;i<N;i++){for(j=0;j<3;j++){a[j][5]=0;for(k=0;k<3;k++){a[j][5]=a[k][4]*a[j][k]+a[j][5];}a[j][5]=a[j][5]+a[j][3];}norm=0;for(k=0;k<3;k++)norm=norm+fabs(a[k][4]-a[k][5]);if(norm<error)break;elsefor(k=0;k<3;k++)a[k][4]=a[k][5];}if(norm<error){printf("计算结果为\n");for(i=0;i<3;i++){printf(" x%d = %.3f\n",i+1,a[i][5]);}}elseprintf("迭代失败\n");}题目1#include "stdio.h"#include "math.h"double f(double x){double ans;ans=exp(x);return(ans);}void main(){double a=-1,b=1,error=0.0001,m=1,h,T0,T,F;int k;h=(b-a)/2;T0=h*(f(a)+f(b));while(1){F=0;for(k=1;k<=pow(2.0,m-1);k++)F=F+f(a+(2*k-1)*h);T=T0/2+h*F;if(fabs(T-T0)<error)break;m++;h=h/2;T0=T;}printf("积分结果为I ≈%f\n",T);}#include "stdio.h"double f(double t,double s){double ans;ans=1+t*t;return(ans);}void main(){double a=0,b=1,h=0.2,x0=0,y0=0,x,k1,k2,k3,y;int N,n;N=(b-a)/h;for(n=1;n<=N;n++){x=x0+h;k1=f(x0,y0);k2=f(x0+h/2,y0+h/2*k1);k3=f(x0+h,y0-h*k1+2*h*k2);y=y0+h/6*(k1+4*k2+k3);printf("第%d次输出结果为(%.8f,%.8f)\n",n,x,y);x0=x;y0=y;}}。
数值分析编程题c语言.word 资料上机实习题一一、题目:已知A与b12.38412,2.115237,-一、题目:已知A与b12.38412,2.115237,:在高斯方法已求出x(m),x(m-1)的基础上,组合新的序列,从而加快收敛速度。
其算式:其中ω是超松弛因子,当ω1时,可以加快收敛速度3 、用消去法求解用追赶消去法求Bx=b的方法:, , , -其算式:其中ω是超松弛因子,当ω1时,可以加快收敛速度3 、用消去法求解用追赶消去法求Bx=b的方法:, , , :y1(0.025000)=0.025000,y2(0.025000)=0.151599,y3(0.025000)=8.33542 6y1(0.045000)=0.045000,y2(0.045000)=0.312860,y3(0.045000)=7.5362 80y1(0.085000)=0.085000,y2(0.085000)=0.560581,y3(0.085000)=4.946 353y1(0.100000)=0.100000,y2(0.100000)=0.628881,y3(0.100000)=4.18 0993四、问题讨论:此程序用高阶方程组的四阶古典的Runge-Kutta方法求解初值问题,控制循环输出语句时,应用时不能用t=0.025(其他判断语句类似),因t是双精度格式的,只能控制t在某一个很小的范围内。
Runge_Kutta方法的优点:1. 精度高,不必用别的方法求开始几点的函数值.2. 可根据f’(t,y)变化的情况与需要的精度自动修改步长.3. 方法稳定.4. 程序简单,存储量少.Runge_Kutta方法的缺点:每步要计算函数值f(t,y)四次,在f(t,y),较复杂时,工作量大, 且每一步缺乏可靠的检查.习题六用对分法求B(即A)的小于23且最接近23的特征值的近似值(误差小于0.1),再用反幂法求B的该特征值的更精确近似值及相应的特征向量。
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 (valuea > 0 && valueb > 0 || valuea <0 && valueb < 0) return;printf("Finding root in the range: [%.3lf, %.3lf]\n", 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:%g\n",x);printf("f(x)=%g\n\n",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("%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 root in the range: [-3.000, -2.000]The x value is:-2.53643f(x)=-0.00124902Finding root in the range: [-1.000, 0.000]The x value is:-0.482857f(x)=0.00012967Finding root in the range: [0.000, 1.000]The x value is:0.482857f(x)=0.00012967Finding root in the range: [2.000, 3.000]The x value is:2.53643f(x)=-0.00124902有时若把判别语句if(fabs(xe)<eps2 || fabs(f(x))<eps1)改为if(fabs(xe)<eps2 && fabs(f(x))<eps1)会提高精度,对同一题运行结果:Finding root in the range: [-3.000, -2.000]The x value is:-2.53644f(x)=-4.26496e-007Finding root in the range: [-1.000, 0.000]The x value is:-0.482861f(x)=-7.3797e-006Finding root in the range: [0.000, 1.000]The x value is:0.482861f(x)=-7.3797e-006Finding root in the range: [2.000, 3.000]The x value is:2.53644f(x)=-4.26496e-007习题2-35. 请用埃特金方法编程求出x=tgx在4.5(弧度)附近的根。
数值分析上机报告姓名:学号:专业:2013年10月27日第一章舍入误差与有效数 设2211NN j S j ==-∑,其精确值为1311221N N ⎛⎫-- ⎪+⎝⎭。
(1)编制按从大到小的顺序22211121311N S N =+++--- ,计算N S 的通用程序。
(2)编制按从小到大的顺序2221111(1)121N S N N =+++---- ,计算NS 的通用程序。
(3)按两种顺序分别计算210S ,410S ,610S ,并指出有效位数。
(编制程序时用单精度) (4)通过本上机题,你明白了什么?解: (1)、(2)题程序见电子版(3)按从大到小顺序:210S =0.740049 有效位数6位 410S =0.749852 有效位数3位 610S =0.749852 有效位数3位 按从小到大顺序:210S =0.740050 有效位数5位 410S =0.749900 有效位数6位610S =0.749999 有效位数6位(4)通过上述实验数据可以看出此次算法使用从小到大的顺序进行得到的数据相对而言更精确,可以得到这样的启示:在计算数值时,要先分析不同算法对结果的影响,避免大数吃小数的现象,找出能得到更精确的结果的算法。
第二章(上机题)Newton 迭代法(1)给定初值0x 及容许误差ε,编制Newton 法解方程()0f x =根的通用程序。
(2)给定方程3()/30f x x x =-=,易知其有三个根1x *=,20x *=,3x *=1.由Newton 方法的局部收敛性可知存在0δ>,当0(,)x δδ∈-时,Newton 迭代序列收敛于根2x *。
试确定尽可能大的δ。
2.试取若干初始值,观察当0(,1)x ∈-∞-,(1,)δ--,(,)δδ-,(,1)δ,(1,)∞时Newton 序列是否收敛以及收敛于哪一个根。
(3)通过本上机题,你明白了什么? 解:(2)1 取esp=0.1,即步长为0.1时,由程序(见电子版)算出max δ=0.78。
本科生课程设计报告实习课程数值分析学院名称管理科学学院专业名称信息与计算科学学生姓名学生学号指导教师实验地点实验成绩二〇一六年六月二〇一六年六摘要常微分方程数值解法是计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题,实用上许多很有价值的常微分方程的解不能用初等函数来表示,常常需要求其数值解.所谓数值解,是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.?关键词:数值解法;常微分1. 实验内容和要求常微分方程初值问题 有精确解2()cos(2)x y x x e x -=+。
要求:分别取步长h = 0.1,0.01,0.001,采用改进的Euler 方法、4阶经典龙格-库塔R -K 方法和4阶Adams 预测-校正方法计算初值问题。
以表格形式列出10个等距节点上的计算值和精确值,并比较他们的计算精确度。
其中多步法需要的初值由经典R-K 法提供。
运算时,取足以表示计算精度的有效位。
2. 算法说明2.1函数及变量说明表1 函数及变量定义1、欧拉方法:1()()(,())i i i i y x y x hf x y x +=+ (i=0,1,2,3,......n-1)(0)y a= (其中a 为初值)2、改进欧拉方法:~1~111()()(,())()()[(,())(,())]2(0)i i i i i i i i i i y x y x hf x y x hy x y x f x y x f x y x y a ++++=+=++=(i=0,1,2......n-1) (其中a 为初值)3、经典K-R 方法: 11213243()6(,)(,)22(,)22(,)i i i i i i i i i i h y y K f x y h hK f x y K h h K f x y K K f x h y hK +⎧=+⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩4、4阶adams 预测-校正方法 预测: 校正:Adsms 内插外插公式联合使用称为Adams 预测-校正系统,利用外插公式计算预测,用内插公式进行校正。
上机实习题一一、题目:已知A 与b12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.0317432.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 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,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.6784392} 1.用Household 变换,把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=1,2,,……,n -2 2 、用超松弛法求解其基本思想:在高斯方法已求出x (m),x (m-1)的基础上,组合新的序列,从而加快收敛速度。
其算式:⎪⎪⎩⎪⎪⎨⎧=+⋅=+-+----=][][0][0][~][]][[][]1[0]][[]1][[]1[]][[]1][[][~i X i X i X i X w i X i i B i b i X i i B i i B i X i i B i i B i X 其中ω是超松弛因子,当ω>1时,可以加快收敛速度 3 、用消去法求解用追赶消去法求Bx=b 的方法:][]1[1i b i d =+ , ]][1[]2[1i i B i a +=+ ,][]1[]1][[][]][[]1[]1][[i b i X i i B i X i i B i X i i B =+⋅++⋅+-⋅-]][[]1[1i i B i b =+ , ]1][[]1[1+=+i i B i c , q1[0]=0 , u1[0]=0 ,8,,2,1]),[1][1][1(][1][1⋅⋅⋅=⋅+-=i i q i a i b i c i q9,,2,1]),1[1][1][1(])[1][1][1(}[1⋅⋅⋅=-⋅+⋅-=i i q i a i b i u i a i d i u x[9]=u1[9]1,,7,8],[1]1[][1][⋅⋅⋅=++⋅=i i u i x i q i x三、1.计算程序:#include "math.h" #include "stdio.h" #define ge 8 void main() {int sign(double x); double a[][9]= {{12.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.101011,-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.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.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 k,h,s,w; int i,j,n,m,g; double u[9],x1[9],y[9],q[9],b1[9][10],x[9]; double b[9]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392}; for(j=0;j<7;++j) /*Household 变换 */ { s=0.0; for(i=j+1;i<9;++i) s=s+a[i][j]*a[i][j]; s=sqrt(s); h=(a[j+1][j]>0)?(s*s+s*a[j+1][j]):(s*s-s*a[j+1][j]); for(g=0;g<9;++g) {if (g<=j)u[g]=0;else if (g==j+1)u[g]=a[j+1][j]+s*sign(a[j+1][j]);else u[g]=a[g][j];}for(m=0;m<9;++m){y[m]=0;for(n=0;n<9;++n)y[m]=y[m]+a[m][n]*u[n];y[m]=y[m]/h;}k=0;for(i=0;i<9;++i)k=k+u[i]*y[i];k=0.5*k/h;for(i=0;i<9;++i)q[i]=y[i]-k*u[i];for(n=0;n<9;++n)for(m=0;m<9;++m)a[m][n]=a[m][n]-(q[m]*u[n]+u[m]*q[n]);}printf("Household:\n");for(i=0;i<9;++i)for(j=0;j<9;++j){if (j%9==0)printf("\n");printf("%-9.5f",a[i][j]);}printf("\n");w=1.4; /*超松弛法*/ for(i=0;i<9;i++)x1[i]=0;for(i=0;i<9;i++)for(j=0;j<9;j++){if(i==j)b1[i][j]=0;else b1[i][j]=-a[i][j]/a[i][i];}for(i=0;i<9;i++)b1[i][9]=b[i]/a[i][i];for(n=0;n<9;n++)for(i=0;i<9;i++){s=0;for(j=0;j<9;j++)s=s+b1[i][j]*x1[j];s=s+b1[i][9];x1[i]=x1[i]*(1-w)+w*s;}for(i=0;i<9;i++){if (i==5)printf("\n");printf("x%d=%-10.6f",i,x1[i]);}printf("\n");u[0]=a[0][0]; /*以下是消去法*/y[0]=b[0];for(i=1;i<9;i++){q[i]=a[i][i-1]/u[i-1];u[i]=a[i][i]-q[i]*a[i-1][i];y[i]=b[i]-q[i]*y[i-1];}x[ge]=y[ge]/u[ge];for(i=ge-1;i>=0;i--)x[i]=(y[i]-a[i][i+1]*x[i+1])/u[i];for(i=0;i<9;i++){if (i==5)printf("\n");printf(" x%d=%-10.6f",i,x[i]);}}int sign(double x){int z;z=(x>=(1e-6)?1:-1);return(z);}2.打印结果:Household:12.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.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.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四、问题讨论:此程序具有很好的通用性。