数值分析程序代码

埃特金#include <iostream.h>#include <math.h>float f(float x){float y;y=exp(-x)return y;}float dd(float x0,float e,float n){float x1,x2;int k;for(k=1;k<=n;k++){x1=f(x0);x2=f(x1);x2=x2-(x2-x1)*(x2-x1)/(x2-2*x1+x0);if(fabs(x2-x0)<e)return x2;cout<<k<<" "<<x2<<endl;x0=x2;}cout<<"error"<<endl;}void main(){float x0,e,n,y;cin>>x0>>e>>n;cout<<endl;y=dd(x0,e,n);cout<<y<<endl;}变步长梯形积分#include <stdio.h>#include <iostream.h>#include <math.h>double a,b,h;double f(double x){double y;y=4/(1+x*x);return y;}double infe(double e,int n){int i;double s,tn,t2n,x;s=(1.0+f(b))/2.0;h=(b-a)/n;x=a+h;for(i=1;i<=n-1;i++){s=s+f(x);x=x+h;}t2n=s*h;dox=a+h/2.0;for(i=1;i<=n;i++){s=s+f(x);x=x+h;}n=2*n;h=h/2.0;t2n=s*h;cout<<t2n<<endl;}while(fabs(t2n-tn)>=e);return t2n;}void main(){ double p,e;int n;cin>>a>>b>>n>>e;cout<<endl;p=infe(e,n);cout<<p<<endl;}变步长新浦生积分#include <stdio.h>#include <iostream.h>#include <math.h>float a,b,e;int n;float f(float x){float y;y=4/(1+x*x);return y;}float simpson(float e,int n) {int i;float s,sn,s2n,p,x,h;h=(b-a)/n;s=f(a)+f(b);p=0;x=a+h;for(i=1;i<=n-1;i++){if((i%2)!=0){p=p+4*f(x);x=x+h;}else{s=s+2*f(x);x=x+h;} }s=s+p;s2n=s*h/3;dox=a+h/2;s=s-p/2;p=0;for(i=1;i<=n;i++){p=p+4*f(x);x=x+h;}s=s+p;n=n*2;h=h/2;s2n=s*h/3;}while(fabs(s2n-sn)>e);cout<<n<<endl;return s2n;}void main(){float t;cin>>a>>b>>e>>n;cout<<endl;t=simpson(e,n);cout<<t<<endl;}迭代#include <iostream.h>#include <math.h>float f(float x){float y;y=exp(-x);return y;}float dd(float x0,float e,float n){float x1;int k;for(k=1;k<=n;k++){x1=f(x0);if(fabs(x1-x0)<e)return x1;cout<<k<<" "<<x1<<endl;x0=x1;}cout<<"error"<<endl;}void main(){float x0,e,n,y;cin>>x0>>e>>n;cout<<endl;y=dd(x0,e,n);cout<<y<<endl;}定步长提醒积分#include <stdio.h>#include <iostream.h>#include <math.h>float fna(float x);void main(){float a,b,s,h;int n,i,j,k,right;float y[100];cin>>a>>b>>n>>right;cout<<endl;h=(b-a)/n;if(right==0){s=(fna(a)+fna(b))/2;for(i=1;i<=n-1;i++)s=s+fna(a+i*h);}else{for(j=0;j<=n;j++)cin>>y[j];s=(y[0]+y[n])/2;for(k=1;k<=n-1;k++)s=s+y[k];}、s=s*h;cout<<s;}float fna(float x){float y;y=x;return y;}二分法#include <iostream.h>float f(float x){float y;y=(x*x*x)-x-1;return y;}float fot(float a,float b,float eps) {float x,y,y0;do{y0=f(a);x=(a+b)/2;y=f(x);if(y*y0>0)a=x;elseb=x;}while(b-a>eps);return x;}void main(){float a,b,eps,x;cout<<"输入范围和误差范围:";cin>>a>>b>>eps;cout<<endl;x=fot(a,b,eps);cout<<x<<endl;}改进欧拉#include <iostream.h>#include <math.h>#include <iomanip.h>double f(double x,double y){double w;w=y-2*x/y;return w;}double fr(double x,double y){double w;w=pow(1+x*2,y);return w;}void main(){double x,y,h,yp,yc,y0;int i,n=10;x=0;y=1;h=0.1;cout<<" x:"<<" y:"<<" y0:"<<endl;for(i=1;i<=n;i++){yp=y+h*f(x,y);yc=y+h*f(x+h,yp);y=(yp+yc)/2;x=x+h;y0=fr(x,0.5);cout<<setw(5)<<setprecision(4)<<x;cout<<setw(10)<<setprecision(7)<<y;cout<<setw(10)<<setprecision(7)<<y0<<endl;} }高斯赛德尔#include <iostream.h>#include <math.h>float a[100][100];float b[100];float x[100];void main(){int n,i,j;cout<<"输入n:";cin>>n;cout<<endl;cout<<"输入a[n][n]"<<endl;for(i=1;i<=n;i++)for(j=1;j<=n;j++)cin>>a[i][j];cout<<"输入b[n]"<<endl;for(i=1;i<=n;i++)cin>>b[i];for(i=1;i<=n;i++){x[i]=0;for(j=1;j<=n;j++)cout<<a[i][j]<<" ";cout<<b[i];cout<<endl;}float t,e,error,q;cout<<"输入e:"<<endl;cin>>e;cout<<endl;do{error=0;for(i=1;i<=n;i++){t=x[i];q=0;for(j=1;j<=n;j++)if(j!=i)q=q+a[i][j]*x[j];x[i]=(b[i]-q)/a[i][i];if(fabs(x[i]-t)>error)error=x[i]-t;}}while(error>e);for(i=1;i<=n;i++)cout<<"x["<<i<<"]="<<x[i]<<endl;}加速迭代#include <stdio.h>#include <iostream.h>#include <math.h>float f(float x);float fr(float x);void main(){float x0,x1,e,q;cin>>x0>>e;x1=x0;do{x0=x1;q=fr(x0);x1=f(x0)+q/(1-q)*(f(x0)-x0);}while(fabs(x1-x0)>=e);cout<<x1<<endl;}float f(float x){float y;x=log(x)+89.62092;y=exp(1.0/3.0*log10(x));return y;}float fr(float x){float y;x=log(x)+89.62092;y=1.0/3.0*exp(-2.0/3.0*log10(x));return y;}科特斯#include <iostream.h>#include <math.h>void main(){double x,a,b,h,s;int k,n;cin>>a>>b>>n;h=(b-a)/n;s=(1-sin(b)/b)*7.0;x=a;for(k=1;k<=n;k++){x=x+h/4.0;s=s+32*sin(x)/x;x=x+h/4.0;s=s+12*sin(x)/x;x=x+h/4.0;s=s+32*sin(x)/x;x=x+h/4.0;s=s+14*sin(x)/x;}s=s*h/90.0;cout<<s<<endl;}拉格朗日插值#include <iostream.h>float x1[100],y1[100];float lagrange(int n,float x){ float y,t;int k,j;y=0;for(k=0;k<=n;k++){ t=1;for(j=0;j<=n;j++)if(j!=k)t=t*((x-x1[j])/(x1[k]-x1[j]));y=y+t*y1[k]; }return y;}void main(){ int i,n;float y,x;cout<<"请输入N:";cin>>n;cout<<endl;cout<<"请输入X1和Y1:";for(i=0;i<=n;i++)cin>>x1[i]>>y1[i];cout<<endl;cout<<"请输入X:";cin>>x;cout<<endl;y=lagrange(n,x);cout<<y<<endl;}龙贝格#include <iostream.h>#include <math.h>float a,b,e;float f(float x){float y;y=4/(1+x*x);return y;}float romberg(float a,float b,float e) {float t1,t2,s1,s2,c1,c2,r1,r2;float x,h,s;int k;h=b-a;t2=(f(a)+f(b))/2*h;k=0;s2=0;c2=0;r2=0;do{t1=t2;s1=s2;c1=c2;r1=r2;k=k+1;x=a+h/2;s=0;while(x<b){s=s+f(x);x=x+h;}h=h/2;t2=t1/2+s*h;s2=t2+(t2-t1)/3;if(k>=2)c2=s2+(s2-s1)/15;if(k>=3)r2=c2+(c2-c1)/63;}while((k<=3)||(fabs(r2-r1)>=e));return r2;}void main(){float w;cin>>a>>b>>e;cout<<endl;w=romberg(a,b,e);cout<<w<<endl;}牛顿切线#include <stdio.h>#include <iostream.h>#include <math.h>float f(float x);float fr(float x);void main(){float x0,x1,e;cin>>x0>>e;x1=x0-f(x0)/fr(x0);while(fabs(x1-x0)>=e){x0=x1;printf("x=%f",x0);x1=x0-f(x0)/fr(x0);printf("=20.7%f",x1);cout<<endl;}}float f(float x){float y;y=x*x-115;return y;}float fr(float x){float y;y=2*x;return y;}欧拉#include <iostream.h>#include <math.h>#include <iomanip.h>double f(double x,double y) {double w;w=y-2*x/y;return w;}double fr(double x,double y){double w;w=pow(1+x*2,y);return w;}void main(){double x,y,h,y0;int i,n=10;x=0;y=1;h=0.1;cout<<" x:"<<" y:"<<" y0:"<<endl;for(i=1;i<=n;i++){y=y+h*f(x,y);x=x+h;y0=fr(x,0.5);cout<<setw(5)<<setprecision(4)<<x;cout<<setw(10)<<setprecision(7)<<y;cout<<setw(10)<<setprecision(7)<<y0<<endl;}}抛物插值#include <iostream.h>void main(){int n,i,t;float x1[100],y1[100],x,y;cout<<"请输入N:";cin>>n;cout<<endl;cout<<"请输入X1和Y1:";for(i=0;i<=n;i++)cin>>x1[i]>>y1[i];cout<<endl;cout<<"请输入X:";cin>>x;cout<<endl;if(x<=x1[1])t=1;if(x>x1[n-1])t=n-1;for(i=1;i<=n-1;i++){if(x1[i-1]<x&&x<=x1[i]){if((x-x1[i-1])>(x1[i]-x))t=i;if((x-x1[i-1])<=(x1[i]-x))t=i-1;}}y=((x-x1[t])*(x-x1[t+1])*(y1[t-1])/(x1[t-1]-x1[t])/(x1[t-1]-x1[t+1]) +(x-x1[t-1])*(x-x1[t+1])*(y1[t])/(x1[t]-x1[t-1])/(x1[t]-x1[t+1])+(x-x1[t-1])*(x-x1[t])*(y1[t+1])/(x1[t+1]-x1[t-1])/(x1[t+1]-x1[t]));cout<<y<<endl;}秦九韶#include <iostream.h>void main(){int n,i;float a[100],x,v;cin>>n;cout<<endl;for(i=1;i<=n;i++){cin>>a[i]; }cout<<endl;for(i=1;i<=n;i++){cout<<a[i];}cout<<endl;cin>>x;cout<<endl;v=a[n];for(i=1;i<n;i++){v=x*v+a[n-i];}cout<<v;}梯形积分递推法#include <iostream.h>#include <math.h>float f(float x){float y;y=sin(x)/x;return y;}void main(){float a,b,e,h,t1,t2,s,x;cin>>a>>b>>e;cout<<endl;h=b-a;t2=(1+f(b))*h/2;do{t1=t2;s=0;x=a+h/2;while(x<b){s=s+f(x);x=x+h;}h=h/2;t2=t1/2+s*h;}while(fabs(t2-t1)>=e);cout<<t2<<endl;}线性插值#include <iostream.h>void main(){int n,i,t;float x1[100],y1[100],x,y;cout<<"请输入N:";cin>>n;cout<<endl;cout<<"请输入X1和Y1:";for(i=0;i<=n;i++)cin>>x1[i]>>y1[i];cout<<endl;cout<<"请输入X:";cin>>x;cout<<endl;if(x<x1[0]){t=1;}if(x>=x1[n]){t=n;}for(i=0;i<=n;i++){if(x1[i]<=x&&x<x1[i+1])t=i+1;}y=y1[t]+(y1[t]-y1[t-1])*(x-x1[t])/(x1[t]-x1[t-1]);cout<<y<<endl;}新浦生#include <iostream.h>#include <math.h>void main(){double x,a,b,h,s;int k,n;cin>>a>>b>>n;h=(b-a)/n;s=1-sin(b)/b;x=a;for(k=1;k<=n;k++){x=x+h/2.0;s=s+4.0*sin(x)/x;x=x+h/2.0;s=s+2.0*sin(x)/x;}s=s*h/6.0;cout<<s<<endl;}。

（一）复化梯形公式例：求121?x dx -=⎰程序：#include "stdio.h"void main(){double a,b,s,h,x;int i,n;a=-1.0;b=1.0;n=10;h=(b-a)/n;x=a;s=x*x/2;for(i=1;i<n;i++){x=x+h;s=s+x*x;}s=s+b*b/2;s=s*h;printf("s=%f\n",s);}结果：s=0.680000（二）复化辛普森公式例：求130?x dx=⎰程序：#include "stdio.h"void main(){double a,b,c,s,h,x,y;int i,n;a=0.0;b=1.0;n=10;s=0.0;h=(b-a)/n;x=a;y=x+h;c=(x+y)/2;for(i=1;i<=n;i++){s=s+x*x*x+4*c*c*c+y*y*y;x=x+h;y=y+h;c=c+h;}s=s*h/6;printf("s=%f\n",s);}结果：s=0.250000（三）复化高斯公式例：求220?x dx=⎰程序：#include <stdio.h>#include <math.h>main(){double a,b,h,s,x1,x2;int i,n;a=0;b=2;n=20;s=0;h=(b-a)/n;for(i=0;i<n;i++){x1=a+i*h+h/2*(1/1.732+1); x2=a+i*h+h/2*(1-1/1.732); s=s+x1*x1*x1+x2*x2*x2; }s=h/2*s;printf("s=%f\n",s);}结果：s=4.000000(四）迭代法例：求x=x2的解。

程序：#include "stdio.h"#include<math.h>main(){double x,xl,y,yl;int i,j;x=0.5;xl=x;y=0.5;yl=y;for(i=0;;i++){x=x*x;if(fabs(xl-x)<0.0001)break;else xl=x;}for(j=0;;j++){y=sqrt(y);if(fabs(yl-y)<0.0001)break;else yl=y;printf("x=%f,y=%f\n",x,y);}结果：x=0.000000，y=0.999915（五）牛顿迭代法y=f(x),求f(x*)=0。

实验一：拉格朗日插值多项式命名（源程序.cpp及工作区.dsw）：lagrange问题：4//Lagrange.cpp#include <stdio.h>#include <conio.h>#define N 4int checkvalid(double x[], int n);void printLag (double x[], double y[], double varx, int n);double Lagrange(double x[], double y[], double varx, int n);void main (){double x[N+1] = {0.4, 0.55, 0.8, 0.9, 1};double y[N+1] = {0.41075, 0.57815, 0.88811, 1.02652, 1.17520};double varx = 0.5;if (checkvalid(x, N) == 1){printf("\n\n插值结果: P(%f)=%f\n", varx, Lagrange(x, y, varx, N));}else{printf("结点必须互异");}getch();}int checkvalid (double x[], int n){int i,j;for (i = 0; i < n; i++){for (j = i + 1; j < n+1; j++){if (x[i] == x[j])//若出现两个相同的结点，返回-1{return -1;}}}return 1;}double Lagrange (double x[], double y[], double varx, int n) {double fenmu;double fenzi;double result = 0;int i,j;printf("Ln(x) =\n");for (i = 0; i < n+1; i++){fenmu = 1;for (j = 0; j < n+1; j++){if (i != j){fenmu = fenmu * (x[i] - x[j]);}}printf("\t%f", y[i] / fenmu);fenzi = 1;for (j = 0; j < n+1; j++){if (i != j){printf("*(x-%f)", x[j]);fenzi = fenzi * (varx - x[j]);}}if (i != n){printf("+\n");}result += y[i] / fenmu * fenzi;}return result;}运行及结果显示：实验二：牛顿插值多项式命名（源程序.cpp及工作区.dsw）：newton_cz问题：4//newton_cz.cpp#include<stdio.h>#include<iostream.h>#include<conio.h>#define N 4int checkvaild(double x[],int n){int i,j;for(i=0;i<n+1;i++){for(j=i+1;j<n+1;j++)if(x[i]==x[j])return -1;}return 1;}void chashang(double x[],double y[],double f[][N+1]){int i,j,t=0;for(i=0;i<N+1;i++){f[i][0]=y[i];//f[0][0]=y[0],f[1][0]=y[1];f[2][0]=y[2];f[3][0]=y[3];f[4][0]=y[4] }for(j=1;j<N+1;j++)// 阶数j{for(i=0;i<N-j+1;i++) //差商个数if[i][j]=(f[i+1][j-1]-f[i][j-1])/(x[t+i+1]-x[i]);//一阶为f[0][1]~f[3][1];二阶为f[0][2]~f[2][2]//三阶为f[0][3]~f[1][3];四阶为f[0][4]t++;}}double compvalue(double t[][N+1],double x[],double varx){int i,j,r=0;double sum=t[0][0],m[N+1]={sum,1,1,1,1};printf("i\tXi\t F(Xi)\t 1阶\t 2阶\t\t3阶\t 4阶\n");printf("--------------------------------------------------------------------------------");for(i=0;i<N+1;i++){r=i;printf("x%d\t%f ",i,x[i]);for(j=0;j<=i;j++){printf("%f ",t[r][j]);r--;}printf("\n");}printf("--------------------------------------------------------------------------------");/**/ printf("N(x) =\n");printf(" %f\n",t[0][0]);for(i=1;i<N+1;i++){for(j=0;j<i;j++)m[i]=m[i]*(varx-x[j]);m[i]=t[0][i]*m[i];sum=sum+m[i];}for(i=1;i<N+1;i++){ printf(" +%f*",t[0][i]);for(j=0;j<i;j++)printf("(x-%f)",x[j]);printf("\n");}return sum;}void main(){double varx,f[N+1][N+1];double x[N+1]={0.4,0.55,0.8,0.9,1};double y[N+1]={0.41075,0.57815,0.88811,1.02652,1.17520};checkvaild(x,N);chashang(x,y,f);varx=0.5;if(checkvaild(x,N)==1){chashang(x,y,f);printf("\n\n牛顿插值结果: P(%f)=%f\n",varx,compvalue(f,x,varx));}elseprintf("输入的插值节点的x值必须互异!");getch();} 运行及结果显示：实验三：自适应梯形公式命名（源程序.cpp 及工作区.dsw ）：autotrap ][n T问题：计算⎰+=10214dx x π的近似值，使得误差（EPS ）不超过610- //autotrap.cpp #include<stdio.h> #include<conio.h> #include<math.h> #define EPS 1e-6 double f(double x) { return 4/(1+x*x); } double AutoTrap(double(*f)(double),double a,double b,double eps) { double t2,t1,sum=0.0; int i,k; t1=(b-a)*(f(a)+f(b))/2; printf("T(%4d)=%f\n",1,t1); for(k=1;k<11;k++) { for(i=1;i<=pow(2,k-1);i++) sum+=f(a+(2*i-1)*(b-a)/pow(2,k)); t2=t1/2+(b-a)*sum/pow(2,k); printf("T(%4d)=%f\n",(int)pow(2,k),t2); if(fabs(t2-t1)<EPS) break; else t1=t2; sum=0.0; } return sum; } void main(){ double s;s=AutoTrap(f,0.0,1.0,EPS);getch();} 运行及结果显示：实验四：龙贝格算法命名（源程序.cpp 及工作区.dsw ）：romberg问题：求⎰+=10214dx x π的近似值，要求误差不超过610-//romberg.cpp #include <stdio.h> #include <conio.h> #include <math.h> #define N 20 #define EPS 1e-6 double f(double x){return 4/(1+x*x);} double Romberg(double a,double b,double (*f)(double),double eps) { double T[N][N],p,h;int k=1,i,j,m=0;T[0][0]=(b-a)/2*(f(a)+f(b));do{p=0;h=(b-a)/pow(2,k-1);for(i=1;i<=pow(2,k-1);i++)p=p+f(a+(2*i-1)*h/2);T[0][k]=T[0][k-1]/2+p*h/2;for(i=1;i<=k;i++){j=k-i;T[i][j]=(pow(4,i)*T[i-1][j+1]-T[i-1][j])/(pow(4,i)-1); }k++; p=fabs(T[k-1][0]-T[k-2][0]);}while(p>=EPS);k--; while(m<=k){for(i=0;i<=m;i++) printf("T(%d,%d)=%f ",i,m-i,T[i][m-i]);m++;printf("\n");}return T[k][0]; } void main() {printf("\n\n 积分结果 = %f",Romberg(0,1,f,EPS)); getch(); } 运行及结果显示：实验五：牛顿切线法问题：求方程01)(3=--=x x x f 在5.1=x ，6.0=x 附近的根（精度31021-⨯=） //newton_qxf.cpp#include <math.h>#include<conio.h>#include <stdio.h>#define MaxK 1000#define EPS 0.5e-3double f(double x){return x*x*x-x-1;}double f1(double x){return 3*x*x-1;}int newton(double (*f)(double), double (*f1)(double), double &x0, double eps) {double xk, xk1;int count = 0;printf("k\txk\n");printf("-----------------------\n");xk = x0;printf("%d\t%f\n", count, xk);do{if((*f1)(xk)==0)return 2;xk1 = xk - (*f)(xk) / (*f1)(xk);if (fabs(xk - xk1) < eps){count++;xk = xk1;printf("%d\t%f\n", count, xk);x0 = xk1;return 1;}count++;xk = xk1;printf("%d\t%f\n", count, xk);}while(count < MaxK);return 0;}void main(){for(int i=0;i<2;i++){double x=0.6;if(i==1)x=1.5;printf("x0初值为%f:\n",x);if (newton(f, f1, x, EPS) == 1){printf("-----------------------\n");printf("the root is x=%f\n\n\n", x);}else{printf("the method is fail!");}}getch();}实验六：牛顿下山法命名（源程序.cpp 及工作区.dsw ）：newton_downhill 问题：求方程01)(3=--=x x x f 在5.1=x ，6.0=x 附近的根（精度31021-⨯=） //newton_downhill.cpp #include <stdio.h> #include <conio.h> #include <math.h> #include <stdlib.h> #define Et 1e-3//下山因子下界 #define E1 1e-3//根的误差限 #define E2 1e-3//残量精度 double f(double x) { return x * x * x - x - 1; } double f1(double x) { return 3 * x * x - 1; } void errormess(int b){char *mess;switch(b){case -1:mess = "f(x)的导数为0!";break;case -2:mess = "下山因子已越界，下山处理失败";break;default:mess = "其他类型错误!";}printf("the method has faild! because %s", mess);}int Newton(double (*f)(double), double (*f1)(double), double &x0) {int k = 0;double t;double xk, xk1;double fxk, fxk_temp;printf("k t xk f(xk)\n");printf("----------------------------------------------------------\n");printf("%-20d", k);xk = x0;printf("%-15f", x0);fxk = (*f)(xk);printf("%-20f", fxk);printf("\n");for (k = 1; ; k++){t = 1;while(1){printf("%-10d", k);printf("%-10f", t);if((*f1)(xk) != 0){xk1 = xk - t * (*f)(xk) / (*f1)(xk);}else{return -1;}printf("%-15f", xk1);fxk_temp = (*f)(xk1);printf("%-20f", fxk_temp);if(fabs(fxk_temp) >fabs(fxk)){t = t / 2;printf("\n");if (t < Et){return -2;}}else{printf("下山成功\n");break;}}if (fabs(xk-xk1)<E1){x0 = xk1;return 1; } xk = xk1; } }void main() {int b;double x0; x0 = 0.6;b = Newton(f, f1, x0); if (b == 1) printf("\nthe root x=%f\n", x0); else errormess(b); getch(); }运行及结果显示：实验七：埃特金加速算法命名（源程序.cpp 及工作区.dsw ）：aitken问题：求方程01)(3=--=x x x f 在5.1=x ，6.0=x 附近的根（精度31021-⨯=） //aitken.cpp#include <math.h> #include <stdio.h> #include <conio.h> #define MaxK 100 #define EPS 0.5e-3double g(double x) {return x * x * x - 1; }int aitken (double (*g)(double), double &x, double eps) {int k;double xk = x, yk, zk, xk1;printf("k xk yk zk xk+1\n");printf("-------------------------------------------------------------------\n"); for (k = 0;k<MaxK; k++) { yk = (*g)(xk); zk = (*g)(yk);xk1 = xk - (yk - xk) * (yk - xk) / (zk - 2 * yk + xk); printf("%-10d%-15f%-15f%-15f%-15f\n", k, xk, yk, zk, xk1); if (fabs(yk-xk)<=eps) { x = xk1; return k + 1; } xk = xk1; }return -1; }void main () {double x = 1.5; int k;k = aitken(g, x, EPS); if (k == -1) printf("迭代次数越界!\n"); else printf("\n 经k=%d 次迭代，所得方程根为:x=%f\n", k, x); getch(); }运行及结果显示：实验八：正割法问题：求方程01)(3=--=x x x f 在5.1=x 附近的根（精度0.5e-8） //ZhengGe.cpp #include <math.h> #include <stdio.h> #include <conio.h>#define MaxK 1000 #define EPS 0.5e-8double f(double x) {return x*x*x-x-1;}int ZhengGe(double (*f)(double), double x0, double &x1, double eps) {printf("k xk f(xk)\n");printf("---------------------------------------------\n");int k;double xk0, xk, xk1;xk0 = x0;printf("%-10d%-15f%-15f\n", 0, x0, (*f)(x0));xk = x1;for (k=1;k<MaxK;k++){if((*f)(xk)-(*f)(xk0)==0)return -1;xk1 = xk - (*f)(xk) / ( (*f)(xk) - (*f)(xk0) ) * (xk - xk0);printf("%-10d%-15f%-15f\n", k, xk, (*f)(xk));if (fabs(xk1 - xk)<=eps){printf("%-10d%-15f%-15f\n\n", k + 1, xk1, (*f)(xk1));printf("---------------------------------------------\n\n");x1 = xk1;return 1;}xk0 = xk;xk = xk1;}return -2;}void main (){double x0 = 1, x1 = 2;if (ZhengGe(f, x0, x1, EPS) == 1){printf("the root is x = %f\n", x1);}else{printf("the method is fail!");}getch();}实验九：高斯列选主元算法命名（源程序.cpp 及工作区.dsw ）：colpivot问题：求解方程组并计算系数矩阵行列式值 ⎪⎩⎪⎨⎧=-+=+-=-+2240532321321321x x x x x x x x x//colpivot.cpp#include <math.h> #include <stdio.h> #include <conio.h> #define N 3static double aa[N][N]={{1,2,-1},{1,-1,5},{4,1,-2}}; static double bb[N+1]={3,0,2};void main() {int i,j;double a[N+1][N+1],b[N+1],x[N+1],det;double gaussl(double a[][N+1],double b[],double x[]); for(i=1;i<=N;i++) { for(j=1;j<=N;j++) a[i][j]=aa[i-1][j-1]; b[i]=bb[i-1]; }det=gaussl(a,b,x); if(det!=0) { printf("\n 方程组的解为:"); for(i=1;i<=N;i++) printf(" x[%d]=%f",i,x[i]); printf("\n\n 系数矩阵的行列式值＝%f",det); }else printf("\n\n 系数矩阵奇异阵，高斯方法失败 !:"); getch(); }double gaussl(double a[][N+1],double b[],double x[])//a传入增广矩阵(有效元素为a[1,1]...a[n,n+1]),x负责传入迭代初值并返回解向量//(有效值为x[1]...x[n]);返回值：系数矩阵行列式值detA{double det=1.0,F,m,temp;int i,j,k,r;void disp(double a[][N+1],double b[]);printf("消元前增广矩阵:\n");disp(a,b);for(k=1;k<N;k++){temp=a[k][k];r=k;for(i=k+1;i<=N;i++){if(fabs(temp)<fabs(a[i][k])){temp=a[i][k];r=i;}//按列选主元，即确定ik}if(a[r][k]==0)return 0;//如果aik,k=0,则A为奇异矩阵，停止计算if(r!=k){for(j=k;j<=N;j++){a[k][j]+=a[r][j];a[r][j]=a[k][j]-a[r][j];a[k][j]-=a[r][j];}b[k]+=b[r];b[r]=b[k]-b[r];b[k]-=b[r];det=-det;//如果ik!=k,则交换[A,b]第ik行与第k行元素printf("交换第%d, %d行:\n",k,r);disp(a,b);}for(i=k+1;i<=N;i++){m=a[i][k]/a[k][k];for(j=1;j<=k;j++)a[i][j]=0.0;for(j=k+1;j<=N;j++)a[i][j]-=m*a[k][j];b[i]-=m*b[k];}printf("第%d次消元:\n",k);disp(a,b);det*=a[k][k];} //大FOR循环结束x[N]=b[N]/a[N][N];for(i=N-1;i>0;i--){F=0.0;for(j=i+1;j<=N;j++)F+=a[i][j]*x[j];x[i]=(b[i]-F)/a[i][i];}det*=a[N][N];return det;}void disp(double a[][N+1],double b[])//显示选主元及消元运算中间增广矩阵结果 {int i,j;for(i=1;i<=N;i++) {for(j=1;j<=N;j++)printf("%10f\t",a[i][j]); printf("%10f\n",b[i]); } }运行及结果显示：实验十：高斯全主元消去算法 命名（源程序.cpp 及工作区.dsw ）：fullpivot问题：利用全主元消去法求解方程组⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----9555.04525.015.0201.0321x x x //fullpivot.cpp#include <math.h> #include <stdio.h> #include <conio.h> #define N 3 double x[N+1];static double aa[N][N]={{0.01,2,-0.5},{-1,-0.5,2},{5,-4,0.5}}; static double bb[N]={-5,5,9};void main(){int i,j;double a[N+1][N+1],b[N+1],x[N+1],det;int t[N+1];//引入列交换保存x向量的各分量的位置顺序double gaussl(double a[][N+1],double b[],double x[],int t[]);for(i=1;i<=N;i++){for(j=1;j<=N;j++) a[i][j]=aa[i-1][j-1];b[i]=bb[i-1];}for(i=1;i<=N;i++)t[i]=i;det=gaussl(a,b,x,t);if(det!=0){ printf("\n方程组的解为:");for(i=N;i>=1;i--){printf(" x[%d]=%f",t[i],x[i]);if(i>1)printf(" -->");}printf("\n\n系数矩阵的行列式值＝%f",det);}else printf("\n\n系数矩阵奇异阵，高斯方法失败!:");getch();}double gaussl(double a[][N+1],double b[],double x[],int t[])//a传入增广矩阵(有效元素为a[1,1]...a[n,n+1]),x负责传入迭代初值并返回解向量//(有效值为x[1]...x[n]);返回值：系数矩阵行列式值detA{double det=1.0,F,m,temp;int i,j,k,r,s;void disp(double a[][N+1],double b[],int x[]);// printf("消元前增广矩阵:\n");//disp(a,b,t);for(k=1;k<N;k++){temp=a[k][k];r=k;s=k;for(i=k;i<=N;i++)for(j=k;j<=N;j++)if(fabs(temp)<fabs(a[i][j])){temp=a[i][j];r=i;s=j;}//选主元，选取ik,jkif(a[r][s]==0)return 0;if(r!=k){for(j=k;j<=N;j++){a[k][j]+=a[r][j];a[r][j]=a[k][j]-a[r][j];a[k][j]-=a[r][j];}b[k]+=b[r];b[r]=b[k]-b[r];b[k]-=b[r];det=-det;printf("交换第%d, %d行:\n",k,r);disp(a,b,t);}if(s!=k){for(i=0;i<=N;i++){a[i][k]+=a[i][s];a[i][s]=a[i][k]-a[i][s];a[i][k]-=a[i][s];}t[k]+=t[s];t[s]=t[k]=t[s];t[k]-=t[s];det=-det;printf("交换第%d, %d列:\n",k,s);disp(a,b,t);}for(i=k+1;i<=N;i++){m=a[i][k]/a[k][k];for(j=1;j<=k;j++)a[i][j]=0.0;for(j=k+1;j<=N;j++)a[i][j]-=m*a[k][j];b[i]-=m*b[k];} //消元计算printf("第%d次消元:\n",k);disp(a,b,t);det*=a[k][k];} //大FOR循环结束x[N]=b[N]/a[N][N];for(i=N-1;i>0;i--){F=0.0;for(j=i+1;j<=N;j++)F+=a[i][j]*x[j];x[i]=(b[i]-F)/a[i][i];}//回代计算det*=a[N][N];return det;}void disp(double a[][N+1],double b[],int x[])//显示选主元及消元运算中间增广矩阵结果{int i,j;for(i=1;i<=N;i++){for(j=1;j<=N;j++)printf("%f\t",a[i][j]);printf("| x%d = %f\n",x[i],b[i]);}}运行及结果显示：实验十一：LU 分解算法命名（源程序.cpp 及工作区.dsw ）：LU问题：利用LU 法求解方程组 ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡472111312613321x x x //LU.cpp#include<math.h> #include<stdio.h> #include<conio.h> #define N 3static aa[N][N]={{3,1,6},{2,1,3},{1,1,1}}; static bb[N]={2,7,4};void main() {int i,j;double a[N+1][N+1],b[N+1];void solve(double a[][N+1],double b[N+1]); int LU(double a[][N+1]); for(i=1;i<=N;i++) { for(j=1;j<=N;j++) a[i][j]=aa[i-1][j-1]; b[i]=bb[i-1]; }if(LU(a)==1){printf("矩阵L如下\n");for(i=1;i<=N;i++){for(j=1;j<=i;j++)if(i==j)printf("%12d ",1);else printf("%12f",a[i][j]);printf("\n");}printf("\n矩阵U如下\n");for(i=1;i<=N;i++){for(j=1;j<=N;j++)if(i<=j)printf("%12f",a[i][j]);else printf("%12s","");printf("\n");}solve(a,b);for(i=1;i<=N;i++)printf("x[%d]=%f ",i,b[i]);printf("\n");}else printf("\nthe LU method failed!\n");getch();}int LU(double a[][N+1])//对N阶矩A阵进行LU分解，结果L、U存放在A的相应位置{int i,j,k,s;double m,n;for(i=2;i<=N;i++)a[i][1]/=a[1][1];for(k=2;k<=N;k++){for(j=k;j<=N;j++){m=0;for(s=1;s<k;s++)m+=a[k][s]*a[s][j];a[k][j]-=m;}if(a[k][k]==0){printf("a[%d][%d]=%d ",k,k,0);return -1;//存在元素akk为0}for(i=k+1;i<=N;i++){n=0;for(s=1;s<k;s++)n+=a[i][s]*a[s][k];a[i][k]=(a[i][k]-n)/a[k][k];}}return 1;//正常结束}void solve(double a[][N+1],double b[N+1])//利用分解的LU求x//回代求解，L和U元素均在A矩阵相应位置；b存放常数列，返回时存放方程组的解{double y[N+1],F;int i,j;y[1]=b[1];for(i=2;i<=N;i++){F=0.0;for(j=1;j<i;j++)F+=a[i][j]*y[j];y[i]=b[i]-F;}b[N]=y[N]/a[N][N];for(i=N-1;i>0;i--){F=0.0;for(j=N;j>i;j--)F+=a[i][j]*b[j];b[i]=(y[i]-F)/a[i][i];}}运行及结果显示：实验十二：Guass-Sediel 算法命名（源程序.cpp 及工作区.dsw ）：GS问题：利用G －S 法求解方程组 ⎪⎩⎪⎨⎧=+--=-+-=--2.453.82102.7210321321321x x x x x x x x x （精度为41021-⨯）//Guass_sediel.cpp#include "iostream.h"#include "math.h"#include "stdio.h"#include<conio.h>#define N 3 //方程的阶数#define MaxK 100 //最大迭代次数#define EPS 0.5e-4 //精度控制static double aa[N][N]={{10,-1,-2},{-1,10,-2},{-1,-1,5}};static double bb[N]={7.2,8.3,4.2};void main(){int i,j;double x[N+1];double a[N+1][N+1],b[N+1];int GaussSeidel(double a[][N+1],double b[],double x[]);for(i=1;i<=N;i++){for(j=1;j<=N;j++)a[i][j]=aa[i-1][j-1];b[i]=bb[i-1];x[i]=0;}if(GaussSeidel(a,b,x)==1){printf("the roots is:");for(i=1;i<=N;i++)printf(" x[%d]=%f ",i,x[i]);printf("\n");}else printf("\nthe G-S method failed!\n");getch();}int GaussSeidel(double a[][N+1],double b[],double x[])//a传入系数矩阵(有效元素为a[1,1]...a[n,n]),b为方程组右边常数列，x传入迭代初值并返回解向量//x**(k+#x)=Gx**(k)+f{int k=1,i,j;double m,max,t[N+1];while(true){printf("k=%d:",k);max=0.0;for(i=1;i<=N;i++){if(a[i][i]==0)return -1;//存在元素akk为0m=0.0;t[i]=x[i];for(j=1;j<=N;j++)if(j!=i)m+=a[i][j]*x[j];x[i]=(b[i]-m)/a[i][i];if(max<fabs(x[i]-t[i]))max=fabs(x[i]-t[i]);printf(" x[%d]=%f ",i,x[i]);}printf("\n");if(max<EPS)return 1;//正常结束elsek++;if(k>MaxK)return -2;//迭代次数越界}}运行及结果显示：实验十三：SOR 算法命名（源程序.cpp 及工作区.dsw ）：sor问题：利用SOR 法求解方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----111141111411114111144321x x x x （精度为51021-⨯） //SOR.cpp#include "math.h"#include "stdio.h"#include "conio.h"#define N 4#define MaxK 1000#define EPS 0.5e-5static double aa[N][N]={{-4,1,1,1},{1,-4,1,1},{1,1,-4,1},{1,1,1,-4}}; static double bb[N]={1,1,1,1};void main(){int i,j;double x[N+1];double a[N+1][N+1],b[N+1];int Sor(double a[][N+1],double b[],double w,double x[]);for(i=1;i<=N;i++){for(j=1;j<=N;j++)a[i][j]=aa[i-1][j-1];b[i]=bb[i-1];x[i]=0;}if(Sor(a,b,1.3,x)==1){printf("the roots is:");for(i=1;i<=N;i++)printf(" x[%d]=%f ",i,x[i]);printf("\n");}else printf("\nthe G-S method failed!\n");getch();}int Sor(double a[][N+1],double b[],double w,double x[]) {int i,j,k;double loft,comb;for(k=0;k<=N;k++)x[k]=0;for(i=1;i<=MaxK && fabs(comb)>EPS;i++){printf("k=%d:\t",i);comb=0;for(k=1;k<=N;k++){loft=b[k];for(j=1;j<N+1;j++)loft-=a[k][j]*x[j];if(a[k][k]==0)return -1;else loft=w*loft/a[k][k];x[k]=x[k]+loft;printf("x[%d]=%10f \t",k,x[k]);if(fabs(loft)>comb)comb=fabs(loft);}printf("\n");}if(i>MaxK)return -2;elsereturn 1;}运行结果（w=1）=。

一、拉格朗日插值#include<stdio.h>#include<stdlib.h>#include<math.h>void Lagrange(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 main(){float x;printf("输入插值点：");scanf("%f",&x);Lagrange(x);}二、牛顿插值#include<stdlib.h>#include<stdio.h>#include<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 main(){float x;printf("输入插值点：");scanf("%f",&x);ND(x);}三、埃尔米特插值#include<stdio.h>#include<stdlib.h>#include<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;}printf("%f\n",H);}void main(){float x;printf("输入插值点：");scanf("%f",&x);Hermite(x);}四、三次样条插值#include <math.h>#include <stdio.h>#include <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++) a[k]=h[k]/(h[k]+h[k-1]);for(k=1;k <N;k++) c[k]=1-a[k];for(k=1;k <N;k++) g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]);c[0]=a[N]=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++)printf( "%9.5f ",s[i]);printf("\n");}五、复合梯形公式#include<stdio.h>#include<stdlib.h>#include<math.h>double FTX(int 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]=sin(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 main(){float a,b;int n;printf("输入区间上，下限：");scanf("%f %f",&a,&b);printf("输入等分区间数：");scanf("%d",&n);FTX(n,a,b);}六、复合辛普森求积公式#include<stdio.h>#include<stdlib.h>#include<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]=sin(x1[i])*x1[i];}for(i=0;i<n;i++){x2[i]=x1[i]+h/2;y2[i]=sin(x2[i])*x2[i];}for(i=1;i<n;i++) f1=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 main(){float a,b;int n;printf("输入区间上，下限：");scanf("%f %f",&a,&b);printf("输入等分区间数：");scanf("%d",&n);FSP(n,a,b);}七、直接三角分解法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){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]=0;L[i][j]=0;}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++){f1=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[0]=b[0];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++)printf("%f ",L[i][j]);printf("\n");}printf("输出U矩阵：\n");for(i=0;i<3;i++){for(j=0;j<3;j++)printf("%f ",U[i][j]);printf("\n");}printf("输出求解结果：\n");for(i=0;i<3;i++)printf("%f ",x[i]);printf("\n");}八、改进的平方法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){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++) printf("%f ",L[i][j]); printf("\n");}printf("输出U矩阵：\n");for(i=0;i<n;i++){for(j=0;j<n;j++) printf("%f ",U[i][j]); printf("\n");}printf("输出求解结果：\n");for(i=0;i<n;i++) printf("%f ",x[i]);printf("\n");}九、追赶法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){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++) printf("%f ",U[i][j]); printf("\n");}printf("输出求解结果：\n");for(i=0;i<n;i++) printf("%f ",x[i]);printf("\n");}十、雅可比迭代法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){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++)printf("%f ",x[i][j]);printf("\n");}printf("\n输出迭代次数：%d\n",k-1);}十一、高斯—塞德尔迭代法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){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++) printf("%f ",x[i][j]); printf("\n");}printf("\n输出迭代次数：%d\n",k-1);}十二、超松弛迭代法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){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;int i,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++) printf("%f ",x[i][j]); printf("\n");}printf("\n输出迭代次数：%d\n",k-1);}数值分析算法程序班级：计算08Q2班姓名：甄彦福。

列主元法function lianzhuyuan(A,b)n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵Ab=zeros(n,1); %矩阵bX=zeros(n,1); %解Xfor i=1:nfor j=1:nA(i,j)=(1/(i+j-1)); %生成hilbert矩阵Aendb(i,1)=sum(A(i,:)); %生成矩阵bendfor i=1:n-1j=i;top=max(abs(A(i:n,j))); %列主元k=j;while abs(A(k,j))~=top %列主元所在行k=k+1;endfor z=1:n %交换主元所在行a1=A(i,z);A(i,z)=A(k,z);A(k,z)=a1;enda2=b(i,1);b(i,1)=b(k,1);b(k,1)=a2;for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵A(s,j)=0;for p=i+1:nA(s,p)=A(s,p)-m*A(i,p);endb(s,1)=b(s,1)-m*b(i,1);endendX(n,1)=b(n,1)/A(n,n); %回代开始for i=n-1:-1:1s=0; %初始化sfor j=i+1:ns=s+A(i,j)*X(j,1);endX(i,1)=(b(i,1)-s)/A(i,i);endX欧拉法clcclear% 欧拉法p=10; %贝塔的取值T=10; %t取值的上限y1=1; %y1的初值r1=1; %y2的初值%输入步长h的值h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0breakendS1=0:T/h;S2=0:T/h;S3=0:T/h;S4=0:T/h;i=1;% 迭代过程for t=0:h:TY=(exp(-t));R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t);y=y1+h*(-y1);y1=y;r=r1+h*(y1-p*r1);r1=r;S1(i)=Y;S2(i)=R;S3(i)=y;S4(i)=r;i=i+1;endt=[0:h:T];% 红线为解析解，'x'为数值解plot(t,S1,'r',t,S3,'x')改进欧拉法clcclearp=10; %贝塔的取值T=10; %t取值的上限y1=1; %y1的初值r1=1; %y2的初值%输入步长h的值h=input('改进欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0breakendS1=0:T/h;S2=0:T/h;S3=0:T/h;S4=0:T/h;i=1;% 迭代过程for t=0:h:TY=(exp(-t));R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t);k1=-y1;l1=y1-p*r1;k2=-(y1+h*k1);l2=y1+h*k1-p*(r1+h*l1);y=y1+h*(k1+k2)/2;r=r1+h*(k1+k2)/2;r1=r;y1=y;S1(i)=Y;S2(i)=R;S3(i)=y;S4(i)=r;i=i+1;endt=[0:h:T];% 红线为解析解，'x'为数值解plot(t,S1,'r',t,S3,'x')高斯-赛德尔function gaosisaideern=input('n='); %阶数tol=input('tol='); %迭代精度A=zeros(n,n);b=zeros(n,1); %生成b向量for i=1:n %给Hilbert矩阵和b向量赋值for j=1:nA(i,j)=(1/(i+j-1));endb(i,1)=sum(A(i,:));endy=zeros(n,1); %迭代解x1=zeros(n,1); %准确解for i=1:ny(i,1)=0; %迭代解赋初值x1(i,1)=1; %生成准确解endk=0;while norm(y-x1)>=tol %精度控制(采用自动步数控制) k=k+1;for i=1:n %迭代开始a1=0;a2=0;for j=1:i-1a1=a1+A(i,j)*y(j,1);endfor j=i+1:na2=a2+A(i,j)*y(j,1);endy(i,1)=(b(i,1)-a1-a2)/A(i,i);endenddisp('迭代步数k')kdisp(y) %显示yend最速下降法function gaosisaideern=input('阶数n='); %阶数tol=input('迭代精度tol='); %迭代精度eps=input('最速下降法eps=');A=zeros(n,n);b=zeros(n,1); %生成b向量for i=1:n %给Hilbert矩阵和b向量赋值for j=1:nA(i,j)=(1/(i+j-1));endb(i,1)=sum(A(i,:));endy=zeros(n,1); %迭代解x1=zeros(n,1); %准确解t=zeros(n,1);r=zeros(n,1);for i=1:ny(i,1)=0; %迭代解赋初值x1(i,1)=1; %生成准确解endr=b-A*y;while norm(r)>=eps; %先进行最速下降法求得进行赛德尔迭代的初始解yt=(r'*r)/(r'*A*r);s1=t*r;y=y+s1;r=b-A*y;endk=0;while norm(y-x1)>=tol %精度控制(采用自动步数控制)k=k+1;for i=1:n %迭代开始a1=0;a2=0;for j=1:i-1a1=a1+A(i,j)*y(j,1);endfor j=i+1:na2=a2+A(i,j)*y(j,1);endy(i,1)=(b(i,1)-a1-a2)/A(i,i);endenddisp('迭代步数k')disp(k)disp(y) %显示y四阶龙格-库塔法clcclearp=10; %贝塔的取值T=10; %t取值的上限y1=1; %y1的初值r1=1; %y2的初值%输入步长h的值h=input('四阶龙格please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0breakendS1=0:T/h;S2=0:T/h;S3=0:T/h;S4=0:T/h;i=1;% 迭代过程for t=0:h:TY=(exp(-t));R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t);k1=-y1;l1=y1-p*r1;k2=-(y1+h*k1/2);l2=y1+h*k1/2-p*(r1+h*l1/2);k3=-(y1+h*k2/2);l3=y1+h*k2/2-p*(r1+h*l2/2);k4=-(y1+h*k3);l4=y1+h*k3-p*(r1+h*l3);y=y1+h*(k1+2*k2+2*k3+k4)/6;r=r1+h*(l1+2*l2+2*l3+l4)/6;r1=r;y1=y;S1(i)=Y;S2(i)=R;S3(i)=y;S4(i)=r;i=i+1;endt=[0:h:T];% 红线为解析解，'x'为数值解plot(t,S1,'r',t,S3,'x')。

拉格朗日插值#include<iostream>#include<cmath>using namespace std;int lagrange(double *px,double *py,double *pcoeff,int iNum) { double *pk,*ptemp_coeff;int i,j,k,kk;pk=new double [iNum];ptemp_coeff=new double [iNum];for(i=0;i<iNum;i++) {pcoeff[i]=0.0;ptemp_coeff[i]=0.0;}for(i=0;i<iNum;i++) {pk[i]=1.0;for(j=0;j<iNum;j++) {if(i!=j) {pk[i]*=px[i]-px[j]; }}}for(i=0;i<iNum;i++) {ptemp_coeff[0]=1.0;for(j=0,k=0;j<iNum;j++) {if(j!=i) {for(k++,kk=k;kk>0;kk--) {ptemp_coeff[kk]-=ptemp_coeff[kk-1]*px[j];}}}for(j=0;j<iNum;j++)pcoeff[j]+=py[i]*ptemp_coeff[j]/pk[i];ptemp_coeff[j]=0.0;}}delete pk,ptemp_coeff;return 0;}double lagcaculate(double *pcoeff,double dxx,int iNum) { double dsum;int i;dsum=0.0;for(i=0;i<iNum;i++) {dsum=dsum*dxx+pcoeff[i];}return dsum;}int main() {int iNum,i,iisrandom;double *px,*py,*pcoeff;cout<<"本程序求解拉格朗日插值求解函数近似值" <<endl<<"请输入已知点个数:"<<endl;cin>>iNum;cout<<"是否进行自动测试（测试函数为正弦函数）,是输入1，否输入0：";cin>>iisrandom;px=new double [iNum];py=new double [iNum];pcoeff=new double [iNum];if(iisrandom==0) {cout<<"请输入iNum个已知点的横坐标:"<<endl;for(i=0;i<iNum;i++) {cin>>px[i];}cout<<"请输入iNum个已知点的纵坐标:"<<endl;for(i=0;i<iNum;i++)cin>>py[i];}else {for(i=0;i<iNum;i++) {px[i]=i;py[i]=sin(i);}}lagrange(px,py,pcoeff, iNum);int testnum;cout<<"输入测试值数量";cin>>testnum;double *pxvalue;pxvalue=new double [testnum];double *pyvalue;pyvalue=new double [testnum];cout<<"输入"<<testnum<<"个需要测试的值"<<endl;for(i=0;i<testnum;i++) {cin>>pxvalue[i];}for(i=0;i<testnum;i++) {pyvalue[i]=lagcaculate(pcoeff,pxvalue[i],iNum);py[i]=sin(pxvalue[i]);}cout<<"拉格朗日公式求得的值"<<" "<<"标准函数值为"<<endl;for(i=0;i<testnum;i++) {cout<<pyvalue[i]<<" "<<py[i]<<endl; }delete px,py,pcoeff,pxvalue,pyvalue;return 0;}牛顿插值#include <iostream.h>#include <math.h>void main() {int n,i,j;double A[50][50];double x[50],y[50];double K=1,X=0,N=0,P;cout<<"请输入所求均差阶数:";cin>>n;for(i=0;i<=n;i++) {cout<<"请输入x"<<i<<"=";cin>>x[i];cout<<"请输入y"<<i<<"=";cin>>y[i];A[i][0]=x[i]; A[i][1]=y[i];}for(j=2;j<=n+1;j++) {for(i=1;i<=n;i++) {A[i][j]=(A[i][j-1]-A[i-1][j-1])/(A[i][0]-A[i-j+1][0]);}}for(i=0;i<=n;i++) {cout<<"输出第"<<i<<"阶均差为："<<A[i][i+1]<<endl;}cout<<"请所要代入计算的x的值：X=";cin>>X;for(i=0;i<n;i++) {K*=X-x[i];N+=A[i+1][i+2]*K;P=A[0][1]+N;}c o u t<<"所要求的函数值为：y="<<P<<endl;}复化梯形公式double f(double x){if(x==0) return 1;else return (sin(x)/x);}double FuhuaTixing(int n,double a,double b){double h = (b-a)/n;double x = a;double s = 0;for(int k=0; k< n-1; k++){x += h;s += f(x);}double T = (f(a)+s*2+f(b))*h/2;return T;}int main(){char ans='n';do{cout<<"请输入积分区间(a,b)："<<endl;double a;double b;cin>>a>>b;cout<<"请输入等分份数n： "<<endl;int n; cin>>n;cout<<"由复化梯形公式球的结果："<<FuhuaTixing(n,a,b)<<endl; cout<<"是否要继续？(y/n)";cin>>ans;}while(ans == 'y');return 0;}改进欧拉#include"stdio.h"#include"iostream"using namespace std;double x,x1,y,y1;int main(){double h;int n;cout<<endl<<"input x0=";cin>>x0;cout<<endl<<"input y0=";cin>>y0;cout<<endl<<"intput h=";cin>>h;cout<<endl<<"intput n=";cin>>n;mend_euler euler(h,n);getch();return 1;}class mend_euler{public:mend_euler(double h,int n);double f(double x,double y);private:double h;int n;};mend_euler::mend_euler(double a,int b){int i=1; h=a; n=b;while(i<=n) {x1=x0+h;y1=y0+h/2*(f(x0,y0)+f(x1,y0+h*f(x0,y0)));cout<<endl;cout<<"x1="<<x1<<" y1="<<y1<<" y="<<x0/(1+x0*x0)<<" e="<<y1-x0/(1+x0*x0) <<endl;i++;x0=x1;y0=y1;}}double mend_euler::f(double x,double y){ return 4*x*y/(1+x*x)+1; }四阶龙格库塔#include <iostream.h> #include <math.h>double f1(double t, double x, double y, double z) { return -x + 6 * y + 2 * z; }double f2(double t, double x, double y, double z) { return -y + 3 * z + 2 * sin(t); }double f3(double t, double x, double y, double z) { return -z + pow(t,2) * exp(-t) + cos(t); } double RungeKutta(double x0, double y0, double z0, double a, double b, int n, double x[], double y[], double z[], double (*f1)(double,double,double,double), double (*f2)(double,double,double,do uble), double (*f3)(double,double,double,double)){double k1, k2, k3, k4double m1, m2, m3, m4double n1, n2, n3, n4double h; double t;int count; count =0;h = (b - a) / n;x[0] = x0;y[0] = y0;z[0] = z0;for(int i = 1;i<=n;i++) {t = a + (i - 1) * h;k1 = f1(t, x[i - 1], y[i - 1], z[i - 1]);m1 = f2(t, x[i - 1], y[i - 1], z[i - 1]);n1 = f3(t, x[i - 1], y[i - 1], z[i - 1]);k2 = f1(t+h/2.0, x[i - 1] + k1 * h / 2.0, y[i - 1] + m1 * h / 2.0, z[i - 1] + n1 * h / 2.0);m2 = f2(t+h/2.0, x[i - 1] + k1 * h / 2.0, y[i - 1] + m1 * h / 2.0, z[i - 1] + n1 * h / 2.0);n2 = f3(t+h/2.0, x[i - 1] + k1 * h / 2.0, y[i - 1] + m1 * h / 2.0, z[i - 1] + n1 * h / 2.0); k3 = f1(t+h/2.0, x[i - 1] + k2 * h / 2.0, y[i - 1] + m2 * h / 2.0, z[i - 1] + n2 * h / 2.0);m3 = f2(t+h/2.0, x[i - 1] + k2 * h / 2.0, y[i - 1] + m2 * h / 2.0, z[i - 1] + n2 * h / 2.0);n3 = f3(t+h/2.0, x[i - 1] + k2 * h / 2.0, y[i - 1] + m2 * h / 2.0, z[i - 1] + n2 * h / 2.0);k4 = f1(t+h/2.0, x[i - 1] + k3 * h / 2.0, y[i - 1] + m3 * h / 2.0, z[i - 1] + n3 * h / 2.0);m4 = f2(t+h/2.0, x[i - 1] + k3 * h / 2.0, y[i - 1] + m3 * h / 2.0, z[i - 1] + n3 * h / 2.0);n4 = f3(t+h/2.0, x[i - 1] + k3 * h / 2.0, y[i - 1] + m3 * h / 2.0, z[i - 1] + n3 * h / 2.0);x[i] = x[i - 1] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;y[i] = y[i - 1] + h * (m1 + 2 * m2 + 2 * m3 + m4) / 6.0;z[i] = z[i - 1] + h * (n1 + 2 * n2 + 2 * n3 + n4) / 6.0; ++count;cout<<"t= "<<t+0.05<<" "<<"x= "<<x[i]<<" "<<"y= "<<y[i]<<" "<<"z= "<<z[i]<<endl;}cout<<"共计算的次数为："<<count<<endl; return 0;}void main() { double x[110]; double y[110];double z[110];RungeKutta(1.0,-1.0,0.0,0.0,5.0,100,x,y,z,f1,f2,f3); }弦截法#include<iostream>#include<math.h>using namespace std;class imaroot{public:virtual float f(float x)=0;virtual float x(float x1,float x2)=0;};class root:public imaroot{public:float f(float x){float y;y=((x-5.0)*x+16.0)*x-80.0;return y;}float x(float x1,float x2){float x;x=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1));return x;}};int main(){root a;float x,c,d;cout<<"输入任意2点横坐标(x1<x2)"<<endl;cin>>c>>d;if(a.f(c)*a.f(d)>=0){do{cout<<"f(x1)与f(x2)同号重新输入"<<endl;cin>>c>>d;}while(a.f(c)*a.f(d)>=0);}if(a.f(c)*a.f(d)<0)x=a.x(c,d);if(a.f(x)*a.f(c)>0){do{x=a.x(x,d);}while(abs(a.f(x))>=0.0001);}else{do{x=a.x(c,x);}while(abs(a.f(x))>=0.0001);}cout<<"方程根是:"<<x<<endl;return 0;}用C语言实现高斯-赛德尔迭代方法，源程序代码：#include "stdlib.h"#include "stdio.h"#include "conio.h"#include "string.h"#include "math.h"#define N 100float Table(int n,float a[N][N],float b[N]){ int i,j; float c[N][N];printf("Please input the matrix A by row!\n"); label:for(i=0;i<n;i++) { printf("Row %d:",i+1); for(j=0;j<n;j++)scanf("%f",&a[i][j]); }printf("Please input the vector b:");for(i=0;i<n;i++)scanf("%f",&b[i]);for(i=0;i<n;i++)for(j=0;j<n;j++) {if(i==j) {c[i][j]=0;continue; }c[i][j]=-a[i][j]/a[i][i]; } printf("\nThe matrix A and vector b:\n"); for(i=0;i<n;i++) { for(j=0;j<n;j++)printf("%10.5f",a[i][j]);printf("%10.5f",b[i]);printf("\n");}printf("\nThe Gauss-Seidel iterative scheme:\n");for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%10.5f",c[i][j]);printf("%10.5f",b[i]/a[i][i]);printf("\n"); } }float init_vec(int n,float x[N]){I nt i;printf("Please input the initial iteration vector x:");for(i=0;i<n;i++)scanf("%f",&x[i]);printf("\nThe initial iteration vector x:\n"); for(i=0;i<n;i++) printf("%10.5f",x[i]);printf("\n");}float gs(int n,float a[N][N],float b[N],float x[N]) {int i,j,k;float tmp1,tmp2,x2[N];for(k=0;k<10001;k++){for(i=0;i<n;i++)x2[i]=x[i];for(i=0;i<n;i++){tmp1=0.0;tmp2=0.0;for(j=0;j<i;j++)tmp1+=a[i][j]*x[j];for(j=i+1;j<n;j++)tmp2+=a[i][j]*x2[j];x[i]=(b[i]-tmp1-tmp2)/a[i][i];}for(i=0,j=0;i<n;i++)if(fabs(x2[i]-x[i])<0.00001) j++;if(j==n){printf("\nThis Gauss-Seidel iterative scheme is convergent!");printf("\nNumber of iterations: %d",k+1);break; }if(k==499){printf("\nThis Jacobi iterative scheme may be not convergent!"); break; }}printf("\nThe results:\n");for(i=0;i<n;i++)printf("%12.7f",x[i]);}int main() {int n;float x[N],a[N][N],b[N];printf("Input n:");scanf("%d",&n);Table(n,a,b);init_vec(n,x);gs(n,a,b,x);getch();}高斯消去#include "Stdio.h"#include "Conio.h"/*L是矩阵的行减1,从程序上看是最外层循环的次数N 对应矩阵的行数，M对应矩阵的列数可以通过改变L、N、M来控制矩的阶数*/#define L 3#define N 4#define M 5void gauss(double a[N][M],double x[N]){int i,j,l,n,m,k=0;double temp[N];/*第一个do-while是将增广矩阵消成上三角形式*/do{n=0;for(l=k;l<L;l++)temp[n++]=a[l+1][k]/a[k][k];for(m=0,i=k;i<N;i++,m++)for(j=k;j<M;j++)a[i+1][j]-=temp[m]*a[k][j];k++;}while(k<N) ;/*第二个do-while是将矩阵消成对角形式，并且重新给k赋值*/k=L-1;do{n=0;for(l=k;l>=0;l--)temp[n++]=a[k-l][k+1]/a[k+1][k+1];for(m=0,i=k;i>=0;i--,m++)for(j=k;j<M;j++)a[k-i][j]-=temp[m]*a[k+1][j];k--;}while(k>=0) ;/*下一个for是解方程组*/for(i=0;i<N;i++)x[i]=a[i][N]/a[i][i];}void menu(){printf("\n _ _ _ _ _\n");printf(" 1.operation\n");printf(" 2.exit");printf("\n _ _ _ _ _\n");}main(){int i,j,choose;double a[N][M]={0},answer[N];clrscr();while(1){leep:menu();scanf("%d",&choose);switch(choose){case 1:printf("!!The size of Maxrix is %d * %d,each line enter %d element:\n ",N,M,M); for(i=0;i<N;i++){printf("Enter the Matrix's %d line:\n",i);for(j=0;j<N+1;j++)scanf("%lf",&a[i][j]);}printf("\nthe corss matrix is:\n_ _ _ _ _\n");gauss(a,answer);for(i=0;i<N;i++){for(j=0;j<M;j++)printf("%-2lf ",a[i][j]);putchar('\n');}printf("_ _ _ _ _\nthe solve is:\n");for(i=0;i<N;i++)printf("x%d=%lf\n",i+1,answer[i]); break;case 2:exit(0);break;default:printf("input error:\n");goto leep;}}getch();}。

function x=tridiagsolver(a,b)[n,n]=size(a);for i=1:nif i==1l(i)=a(i,i);y(i)=b(i)/l(i);elsel(i)=a(i,i)-a(i,i-1)*u(i-1);y(i)=(b(i)-y(i-1)*a(i,i-1))/l(i);endif i<nu(i)=a(i,i+1)/l(i);endendx(n)=y(n);for j=n-1:-1:1x(j)=y(j)-u(j)*x(j+1);end拉格朗日function yh=lagrange(x,y,xh)n=length(x);m=length(xh);yh=zeros(1,m);c1=ones(n-1,1);c2=ones(1,m);for i=1:nxp=x([1:i-1 i+1:n]);yh=yh+y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));end线性x=[x1,x2] y=[y1.y2] xh=[xh]抛物线X=[x1,x2,x3] y=[y1，y2，y3] xh=[xh]牛顿差商（输入x,y为列向量）function [p,q]=chashang(x,y)n=length(x);p(:,1)=x;p(:,2)=y;for j=3:n+1p(1:n+2-j,j)=diff(p(1:n+3-j,j-1))./(x(j-1:n)-x(1:n+2-j)); endq=p(1,2:n+1)';三次样条x=[0 1 2 3];y=[0.2 0 0.5 2.0 1.5 -1];pp=csape(x,y,'complete')[breaks,coefs,npolys,ncoefs,dim]=unmkpp(pp)x=[0.24 0.65 0.95 1.24 1.73 2.01 2.23 2.52 2.77 2.99]';y=[0.23 -0.26 -1.1 -0.45 0.27 0.1 -0.29 0.24 0.56 1]';A=[log(x) cos(x) exp(x)];Z=A\y;a0=Z(1)a1=Z(2)a2=Z(3)x=[0 0.25 0.50 0.75 1.00];y=[1.00 1.284 1.6487 2.1170 2.7183];p=polyfit(x,y,2)a2=p(1)a1=p(2)a0=p(3)复合中点function I=fmid(fun,a,b,n)h=(b-a)/n;x=linspace(a+h/2,b-h/2,n);y=feval(fun,x);I=h*sum(y);复合梯形function I=ftrapz(fun,a,b,n)h=(b-a)/n;x=linspace(a,b,n+1);y=feval(fun,x);I=h*(0.5*y(1)+0.5*y(n+1)+sum(y(2:n)));复合辛普森function I=fsimpson(fun,a,b,n)h=(b-a)/n;x=linspace(a,b,2*n+1);y=feval(fun,x);I=h/6*(y(1)+y(2*n+1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n)));雅克比迭代function [x,iter]=jacobi(A,b,tol)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x=zeros(size(b));for iter=1:500x=D\(b+U*x+L*x);error=norm(b-A*x)/norm(b);if (error<tol)break;endGS迭代function [x,iter]=GS(A,b,tol)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x=zeros(size(b));for iter=1:500x=(D-L)\(b+U*x);error=norm(b-A*x)/norm(b);if (error<tol)break;endendSOR迭代function [x,iter]=SOR(A,b,omega,tol)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x=zeros(size(b));for iter=1:500x=(D-omega*L)\(omega*b+(1-omega)*D*x+omega*U*x); error=norm(b-A*x)/norm(b);if (error<tol)break;endend二分法function v=f(x)v=x^3-x-1;function x=erfenfa(f,a,b,tol)if nargin<4tol=1e-5;endfa=feval(f,a);fb=feval(f,b);while abs(a-b)>tolx=(a+b)/2;fx=feval(f,x);if sign(fx)==sign(fa)a=x;fa=fx;elseif sign(fx)==sign(fb)b=x;fb=fx;endend>> x=erfenfa('f',1,2)牛顿法function [x,it]=newton(f,g,x0,tol)it=0;done=0;while ~donex=x0-feval(f,x0)/feval(g,x0);it=it+1;done=(norm(x-x0)<=tol);if ~donex0=x;endendfunction r=f(x)r=polyval([1,2,10,-20],x);function r=g(x)p=polyder([1,2,10,-20]);r=polyval(p,x);牛顿下山法乘幂法function [t,y]=chengmifa(a,xinit,ep) v0=xinit;[tv,ti]=max(abs(v0));lam0=v0(ti);u0=v0/lam0;flag=0;while ~flagv1=a*u0;[tv,ti]=max(abs(v1));lam1=v1(ti);u0=v1/lam1;err=abs(lam0-lam1);if err<epflag=1endlam0=lam1;endt=lam1;y=u0;反幂法function [t,y]=fanmifa(a,xinit,ep)v0=xinit;[tv,ti]=max(abs(v0));lam0=v0(ti);u0=v0/lam0;flag=0;while ~flagv1=inv(a)*u0;[tv,ti]=max(abs(v1));lam1=v1(ti);u0=v1/lam1;err=abs(lam0-lam1);if err<epflag=1;endlam0=lam1;endt=1/lam1;y=u0;改进欧拉法function [x,y]=odeIEuler(f,y0,a,b,n) h=(b-a)/n;x=a:h:b;y(1)=y0;for i=1:nyp=y(i)+h*feval(f,x(i),y(i));yc=y(i)+h*feval(f,x(i+1),yp);y(i+1)=1/2*(yp+yc);endfunction r=f(x,y)r=y-2*x/y;。

1、%用牛顿法求f(x)=x-sin x 的零点，e=10^(-6)disp('牛顿法');i=1;n0=180;p0=pi/3;tol=10^(-6);for i=1:n0p=p0-(p0-sin(p0))/(1-cos(p0));if abs(p-p0)<=10^(-6)disp('用牛顿法求得方程的根为')disp(p);disp('迭代次数为：')disp(i)break;endp0=p;endif i==n0&&~(abs(p-p0)<=10^(-6))disp(n0)disp('次牛顿迭代后无法求出方程的解')end2、disp('Steffensen加速');p0=pi/3;for i=1:n0p1=0.5*p0+0.5*cos(p0);p2=0.5*p1+0.5*cos(p1);p=p0-((p1-p0).^2)./(p2-2.*p1+p0);if abs(p-p0)<=10^(-6)disp('用Steffensen加速求得方程的根为')disp(p);disp('迭代次数为：')disp(i)break;endp0=p;endif i==n0&&~(abs(p-p0)<=10^(-6))disp(n0)disp('次Steffensen加速后无法求出方程的解')end1、%使用二分法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，%误差限为 e=10^-4disp('二分法')a=0.2;b=0.26;tol=0.0001;n0=10;fa=600*(a.^4)-550*(a.^3)+200*(a.^2)-20*a-1;for i=1:n0p=(a+b)/2;fp=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1;if fp==0||(abs((b-a)/2)<tol)disp('用二分法求得方程的根p=')disp(p)disp('二分迭代次数为：')disp(i)break;endif fa*fp>0a=p;else b=p;endendif i==n0&&~(fp==0||(abs((b-a)/2)<tol))disp(n0)disp('次二分迭代后没有求出方程的根')end2、%使用牛顿法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，%误差限为 e=10^-4disp('牛顿法')p0=0.3;for i=1:n0p=p0-(600*(p0.^4)-550*(p0.^3)+200*(p0.^2)-20*p0-1)./(2400*(p0.^3) -1650*p0.^2+400*p0-20);if(abs(p-p0)<tol)disp('用牛顿法求得方程的根p=')disp(p)disp('牛顿迭代次数为：')disp(i)break;endp0=p;endif i==n0&&~(abs(p-p0)<tol)disp(n0)disp('次牛顿迭代后没有求出方程的根')end3、%使用割线法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，%误差限为 e=10^-4disp('割线法')p0=0.2;p1=0.25;q0=600*(p0.^4)-550*(p0.^3)+200*(p0.^2)-20*p0-1;q1=600*(p1.^4)-550*(p1.^3)+200*(p1.^2)-20*p1-1;for i=2:n0p=p1-q1*(p1-p0)/(q1-q0);if abs(p-p1)<toldisp('用割线法求得方程的根p=')disp(p)disp('割线法迭代次数为：')disp(i)break;endp0=p1;q0=q1;pp=p1;p1=p;q1=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1;endif i==n0&&~(abs(p-pp)<tol)disp(n0)disp('次割线法迭代后没有求出方程的根')end4、%使用试位法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，%误差限为 e=10^-4disp('试位法')p0=0.2;p1=0.25;q0=600*(p0.^4)-550*(p0.^3)+200*(p0.^2)-20*p0-1;q1=600*(p1.^4)-550*(p1.^3)+200*(p1.^2)-20*p1-1;for i=2:n0p=p1-q1*(p1-p0)/(q1-q0);if abs(p-p1)<toldisp('用试位法求得方程的根p=')disp(p)disp('试位法迭代次数为：')disp(i)break;endq=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1;if q*q1<0p0=p1;q0=q1;endpp=p1;p1=p;q1=q;endif i==n0&&~(abs(p-pp)<tol)disp(n0)disp('次试位法迭代后没有求出方程的根')end5、%使用muller方法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，%误差限为 e=10^-4disp('muller法')x0=0.1;x1=0.2;x2=0.25;h1=x1-x0;h2=x2-x1;d1=((600*(x1.^4)-550*(x1.^3)+200*(x1.^2)-20*x1-1)-(600*(x0.^4)-55 0*(x0.^3)+200*(x0.^2)-20*x0-1))/h1;d2=((600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)-(600*(x1.^4)-55 0*(x1.^3)+200*(x1.^2)-20*x1-1))/h2;d=(d2-d1)/(h2+h1);for i=3:n0b=d2+h2*d;D=(b*b-4*(600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)*d)^0.5;if(abs(d-D)<abs(d+D))E=b+D;else E=b-D;endh=-2*(600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)/E;p=x2+h;if abs(h)<toldisp('用muller方法求得方程的根p=')disp(p)disp('muller方法迭代次数为：')disp(i)break;endx0=x1;x1=x2;x2=p;h1=x1-x0;h2=x2-x1;d1=((600*(x1.^4)-550*(x1.^3)+200*(x1.^2)-20*x1-1)-(600*(x0.^4)-55 0*(x0.^3)+200*(x0.^2)-20*x0-1))/h1;d2=((600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)-(600*(x1.^4)-55 0*(x1.^3)+200*(x1.^2)-20*x1-1))/h2;d=(d2-d1)/(h2+h1);endif i==n0%条件有待商榷？！disp(n0)disp('次muller方法迭代后没有求出方程的根')end1、%观察Lagrange插值的Runge现象x=-1:0.05:1;y=1./(1+25.*x.*x);plot(x,y),grid on;n=5;x=-1:2/n:1;y=1./(1+25.*x.*x);for i=1:n+1q(1,i)=y(i);endh=0.05;z=-1:h:1;for k=1:2/h+1for i=2:n+1for j=2:iq(j,i)=((z(k)-x(i-j+1))*q(j-1,i)-(z(k)-x(i))*q(j-1,i-1))/(x(i)-x( i-j+1));endendw(k)=q(n+1,n+1);endhold on, plot(z,w,'r'),grid on;%**** n=10 ****n=10;x=-1:2/n:1;y=1./(1+25.*x.*x);for i=1:n+1q(1,i)=y(i);endh=0.05;z=-1:h:1;for k=1:2/h+1for i=2:n+1for j=2:iq(j,i)=((z(k)-x(i-j+1))*q(j-1,i)-(z(k)-x(i))*q(j-1,i-1))/(x(i)-x( i-j+1));endendw(k)=q(n+1,n+1);endhold on,plot(z,w,'k'),grid on;legend ('原始图','n=5','n=10');2、%固支样条插植%********第一段********x=[1,2,5,6,7,8,10,13,17];a=[3,3.7,3.9,4.2,5.7,6.6,7.1,6.7,4.5];n=numel(a);for i=1:n-1h(i)=x(i+1)-x(i);endA=[2*h(1),h(1),0,0,0,0,0,0,0;h(1),2*(h(1)+h(2)),h(2),0,0,0,0,0,0;0,h(2),2*(h(2)+h(3)),h(3),0,0,0,0,0;0,0,h(3),2*(h(3)+h(4)),h(4),0,0,0,0;0,0,0,h(4),2*(h(4)+h(5)),h(5),0,0,0;0,0,0,0,h(5),2*(h(5)+h(6)),h(6),0,0;0,0,0,0,0,h(6),2*(h(6)+h(7)),h(7),0;0,0,0,0,0,0,h(7),2*(h(7)+h(8)),h(8);0,0,0,0,0,0,0,h(8),2*h(8)];e=[3*(a(2)-a(1))/h(1)-3;3*(a(3)-a(2))/h(2)-3*(a(2)-a(1))/h(1);3*(a(4)-a(3))/h(3)-3*(a(3)-a(2))/h(2);3*(a(5)-a(4))/h(4)-3*(a(4)-a(3))/h(3);3*(a(6)-a(5))/h(5)-3*(a(5)-a(4))/h(4);3*(a(7)-a(6))/h(6)-3*(a(6)-a(5))/h(5);3*(a(8)-a(7))/h(7)-3*(a(7)-a(6))/h(6);3*(a(9)-a(8))/h(8)-3*(a(8)-a(7))/h(7);3*(-0.67)-3*(a(9)-a(8))/h(8)];c=inv(A)*e;for i=1:8b(i)=(a(i+1)-a(i))/h(i)-h(i)*(2*c(i)+c(i+1))/3;d(i)=(c(i+1)-c(i))/(3*h(i));endfor i=1:8z=x(i):0.05:x(i+1);w=a(i)+b(i).*(z-x(i))+c(i).*(z-x(i)).^2+d(i).*(z-x(i)).^3; grid on, plot(z,w),hold on;end%********第二段********x=[17,20,23,24,25,27,27.7];a=[4.5,7,6.1,5.6,5.8,5.2,4.1];for i=1:6h(i)=x(i+1)-x(i);endA=[2*h(1),h(1),0,0,0,0,0;h(1),2*(h(1)+h(2)),h(2),0,0,0,0;0,h(2),2*(h(2)+h(3)),h(3),0,0,0;0,0,h(3),2*(h(3)+h(4)),h(4),0,0;0,0,0,h(4),2*(h(4)+h(5)),h(5),0;0,0,0,0,h(5),2*(h(5)+h(6)),h(6)0,0,0,0,0,h(6),2*h(6)];e=[3*(a(2)-a(1))/h(1)-3*3;3*(a(3)-a(2))/h(2)-3*(a(2)-a(1))/h(1);3*(a(4)-a(3))/h(3)-3*(a(3)-a(2))/h(2);3*(a(5)-a(4))/h(4)-3*(a(4)-a(3))/h(3);3*(a(6)-a(5))/h(5)-3*(a(5)-a(4))/h(4);3*(a(7)-a(6))/h(6)-3*(a(6)-a(5))/h(5);3*(-4)-3*(a(7)-a(6))/h(6)];c=inv(A)*e;for i=1:6b(i)=(a(i+1)-a(i))/h(i)-h(i)*(2*c(i)+c(i+1))/3;d(i)=(c(i+1)-c(i))/(3*h(i));endfor i=1:6z=x(i):0.05:x(i+1);w=a(i)+b(i).*(z-x(i))+c(i).*(z-x(i)).^2+d(i).*(z-x(i)).^3; grid on, plot(z,w),hold on;end%********第三段********x=[27.7,28,29,30];a=[4.1,4.3,4.1,3];for i=1:3h(i)=x(i+1)-x(i);endA=[2*h(1),h(1),0,0;h(1),2*(h(1)+h(2)),h(2),0;0,h(2),2*(h(2)+h(3)),h(3);0,0,h(3),2*h(3)];e=[3*(a(2)-a(1))/h(1)-3*0.33;3*(a(3)-a(2))/h(2)-3*(a(2)-a(1))/h(1);3*(a(4)-a(3))/h(3)-3*(a(3)-a(2))/h(2);3*(-1.5)-3*(a(4)-a(3))/h(3)];c=inv(A)*e;for i=1:3b(i)=(a(i+1)-a(i))/h(i)-h(i)*(2*c(i)+c(i+1))/3;d(i)=(c(i+1)-c(i))/(3*h(i));endfor i=1:3z=x(i):0.05:x(i+1);w=a(i)+b(i).*(z-x(i))+c(i).*(z-x(i)).^2+d(i).*(z-x(i)).^3; grid on, plot(z,w),hold on;endgrid on,title('注：横纵坐标的比例不一样！！！');1、%用不动点迭代法求方程 x-e^x+4=0的正根与负根，误差限是10^-6%disp('不动点迭代法');n0=100;p0=-5;for i=1:n0p=exp(p0)-4;if abs(p-p0)<=10^(-6)if p<0disp('|p-p0|=')disp(abs(p-p0))disp('不动点迭代法求得方程的负根为:')disp(p);break;elsedisp('不动点迭代法无法求出方程的负根.')endelsep0=p;endendif i==n0disp(n0)disp('次不动点迭代后无法求出方程的负根')endp1=1.7;for i=1:n0pp=exp(p1)-4;if abs(pp-p1)<=10^(-6)if pp>0disp('|p-p1|=')disp(abs(pp-p1))disp('用不动点迭代法求得方程的正根为')disp(pp);elsedisp('用不动点迭代法无法求出方程的正根');endbreak;elsep1=pp;endendif i==n0disp(n0)disp('次不动点迭代后无法求出方程的正根')end2、%用牛顿法求方程 x-e^x+4=0的正根与负根，误差限是10^-6 disp('牛顿法')n0=80;p0=1;for i=1:n0p=p0-(p0-exp(p0)+4)/(1-exp(p0));if abs(p-p0)<=10^(-6)disp('|p-p0|=')disp(abs(p-p0))disp('用牛顿法求得方程的正根为')disp(p);break;elsep0=p;endendif i==n0disp(n0)disp('次牛顿迭代后无法求出方程的解')endp1=-3;for i=1:n0p=p1-(p1-exp(p1)+4)/(1-exp(p1));if abs(p-p1)<=10^(-6)disp('|p-p1|=')disp(abs(p-p1))disp('用牛顿法求得方程的负根为')disp(p);break;elsep1=p;endendif i==n0disp(n0)disp('次牛顿迭代后无法求出方程的解')end1、使用欧拉法、改进欧拉法和四阶R-K方法求下列微分方程的解。

一、非线性方程求解二分法# include "math.h"# define f(x) sin(x)-x*x/2main(){ float a,b,c,e;int k=0;printf("\n\n\n 二分法求非线性方程的解"); Use Dichotomy solving nonlinear equations printf("\n\n f(x)已在源程序中定义!"); f (x) is defined in the sourceprintf("\n 已知f(x)在[a,b]上连续"); Known f (x) in [a, b] is continuousprintf("\n 请输入a:");scanf("%f",&a);printf(" 请输入b:");scanf("%f",&b);printf(" 请输入允许误差e:");scanf("%f",&e);loop:c=(a+b)/2;printf("\n\n a%d=%f",k,a);printf("\n b%d=%f",k,b);printf("\n x%d=%f",k,c);k++;if((f(a))*(f(c))>0) a=c;else b=c;if (b-a<=e){ printf("\n 最终结果是:");printf("\n x%d=%f",k,(a+b)/2);exit(0);}else goto loop;}运行结果：1.x18=1.4044132.x17=1.324717牛顿法# include "math.h"# define f(x) x*exp(x)-1# define f1(x) exp(x)+x*exp(x)main(){ float a, x0,x1,e1,e2;int N,k;printf("\n\n\n 用牛顿迭代法求非线性方程的解");printf("\n\n\n函数f(x)及f'(x)已在源程序中定义!");printf("\n\n 请输入初值x0:");scanf("%f",&x0);printf(" 输入误差e1:");scanf("%f",&e1);printf(" 输入误差e2:");scanf("%f",&e2);printf(" 输入最大迭代次数N:");scanf("%d",&N);k=1;loop:if(f(x0)==0) {printf("\n %f为所求得根之一!",x0);exit(0);}else if (fabs(f(x0))<e1) { printf("\n奇异标志!建议重取初值或取更小的e1,e2!"); exit(0);} else x1=x0-(f(x0))/(f1(x0));printf("\n 第%d次迭代:",k);printf("\n x%d=%f",k,x1);if (fabs(x1-x0)<e2) { printf("\n %f 为所求近似解!",x1);exit(0);}else if(k==N) { printf("\n迭代失败!建议取更大的N值或另选初值!"); exit(0); }else { x0=x1;k++;}goto loop ;}运行结果：1)xe x－1＝0 x0＝0.5x3＝x4＝0.5671432)x3－x－1＝0 x0＝1x4=x5＝1.3247183)（x-1）2(2x-1)=0x0=0.45 x3=x4＝0.500000x0=0.65 x8=0.5000004 x4=0.500000二、Gauss列主元消去法# define NUMBER 20# include "math.h"float A[NUMBER][NUMBER+1] ,ark;int flag,n;exchange(int r,int k){ int i;for(i=1;i<=n+1;i++)A[0][i]=A[r][i];for(i=1;i<=n+1;i++)A[r][i]=A[k][i];for(i=1;i<=n+1;i++)A[k][i]=A[0][i];}float max(int k){ int i;float temp=0;for(i=k;i<=n;i++)if(fabs(A[i][k])>temp){ temp=fabs(A[i][k]);flag=i;}return temp;}main(){ float x[NUMBER],b[NUMBER];int r,k,i,j;re: printf("\n\n\n 用Gauss列主元消元法解线性方程组");printf("\n input n=");scanf("%d",&n);printf("\n\n 现在输入系数矩阵A和向量b:\n\n");for(i=1;i<=n;i++)for(j=1;j<=n+1;j++){ if(j==n+1) printf(" 请输入b%d:",i);else printf(" 请输入a%d%d:",i,j);scanf("%f",&A[i][j]);}k=1;loop:ark=max(k);if(ark==0){ printf("\n不是合法的方程组!");exit(0);}else if(flag!=k)exchange(flag,k);for(i=k+1;i<=n;i++)for(j=k+1;j<=n+1;j++)A[i][j]=A[i][j]-A[i][k]*A[k][j]/A[k][k];if(k<n-1){ k++; goto loop;}else{ x[n]=A[n][n+1]/A[n][n];for( k=n-1;k>=1;k--){float me=0;for(j=k+1;j<=n;j++){ me=me+A[k][j]*x[j];}x[k]=(A[k][n+1]-me)/A[k][k];}}for(i=1;i<=n;i++)printf("\n x%d=%f",i,x[i]);goto re;}运行结果：1.x1= -0.491058x2= -0.050886x3= 0.3672572.x1=2.946429x2=0.607143x3=-0.142857若不选主元，则1的结果溢出，2不变，原因是1中a11太小，在消元过程中作了分母。