数值分析编程题c语言
数值分析作业(C++)

数值分析上机题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的值,本题为了编程方便,已经把数据写入程序中。
东南大学数值分析上机题C语言编程源程序(前三章)

第一章作业#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]);}。
数值分析编程题c语言

上机实习题一一、题目:已知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程序

b[k*n+i] = t;
}
}
// 为了使a的第k行的第k元素为1,a和b的第k行都除以a(k,k)
t = 1.0 / a[k*n+k];
for(i=k; i<n; i++) a[k*n+i] *= t;
for(i=0; i<n; i++) b[k*n+i] *= t;
for(i=0; i<n; i++)
{
if(i != k && a[i*n+k] != 0.0)
{
t = a[i*n+k];
for(j=k+1; j<n; j++) a[i*n+j] -= a[k*n+j] * t;
for(j=0; j<n; j++) b[i*n+j] -= b[k*n+j] * t;
// 带有列主元的约当消元法
// 功能: 求解线性方程组 Ax = b
// 参数: A - 指向n*n系数矩阵的指针
// b - 常数向量的指针
// n - 方程组的维数
// 返回值:0 - 如果成功。线性方程组的解保存在 b 中
// 1 - 求解失败
// ----------------------------------------------------------------------------
if(j != k) // 交换方程的第 j 行和 k 行
{
for(i=k; i<n; i++)
数值分析五个题目的C语言及MATLAB程序

1.7917595,1.9459101,2.079445,2.1972246,2.3025851}; //函数声明 void canshu(double *h,double *r,double *u,double *d); void wanju(double *u,double *r,double *d,double *m); double yths(double *m,double *h,double p); double ytdhs(double *m,double *h,double p); //主函数 void main() {
break; end end 实验二 C 程序 #include"stdio.h" #include"math.h" //已知量 double x[10]={1,2,3,4,5,6,7,8,9,10}; double fd1=1,fd2=0.1; double fx[10]={0,0.69314718,1.0986123,1.3862944,1.6094378,
} Matlab 程序 function[h,r,u,d]=canshu(x,fx,h,r,u,d,fd1,fd2) %UNTITLED10 Summary of this function goes here % Detailed explanation goes here for i=1:9
h(i)=x(i+1)-x(i); a(i)=(fx(i)-fx(i+1))/(x(i)-x(i+1)); j=i-1;
end end function [q] =ytdhs(m,h,p,fx,x) %UNTITLED14 Summary of this function goes here % Detailed explanation goes here for i=1:10
C语言编程题带答案

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 ,用于存储数组元素的总和。
北航数值分析-实习作业1(C语言详细注释)
《数值分析》计算实习作业《一》北航第一题 设有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,再重复上述步骤即可求得两个特征值。
数值分析算法C语言程序
拉格朗日插值#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 班姓名:甄彦福。
(完整word版)数值分析作业(C语言编程实现)(word文档良心出品)
#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语言
Y=Y+r*y[i];}
printf("%lf\n",Y);
}
2牛顿插值:
#include<stdio.h>
#define N 6
void main()
{
int i,k,n=N-1;
float X,Y,x[N],y[N],c[N],b[N];
printf("请输入xi:\n");
for(i=0;i<N;i++)
int k,N;
double x0,x1;
printf("请输入x0,N\n");
scanf("%lf,%d",&x0,&N);
for(k=1;k<=N;k++)
{
if(fd(x0)==0){printf("奇异标志\n");break;}
x1=x0-(f(x0)/fd(x0));
printf("X%d= %lf\n",k,x1);
}
}
float f(float x,float y)
{
float s;
s=y-(2*x)/y;
return s;
}
5牛顿迭代法:
#include<stdio.h>
#include<math.h>
#define e 1e-7
void main()
{
double f(double x);
double fd(double x);
1拉格朗日插值:
#include<stdio.h>
#define N 3
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
上机实习题一一、题目:已知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四、问题讨论:此程序具有很好的通用性。