数值分析程序代码
数值分析实验代码
埃特金#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）=。
数值分析 算法C语言程序
一、拉格朗日插值#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班姓名：甄彦福。
数值分析的MATLAB程序
列主元法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;。
数值分析实验程序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
数值分析matlab代码
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太小，在消元过程中作了分母。
数值分析-代码
注：凡“in.open("C:\\Users\\Administrator\\Desktop\\data-拉格朗日插值.txt")”等里面的路径自己合理修改。
1-1.拉格朗日插值（用一系列已知点求近似函数，并估计未知点的函数值）#include<fstream.h>#include<stdlib.h>int n; //插值节点个数double *X,*Y;//插值节点数组double x;//估计点void read(){double temp;ifstream in;in.open("C:\\Users\\Administrator\\Desktop\\data-拉格朗日插值.txt");if(!in){cout<<"can not open the file !"<<endl;exit(0);}in>>n;X=new double[n];Y=new double[n];for(int i=0;i<n;i++){in>>temp;X[i]=temp;cout<<temp<<" ";in>>temp;Y[i]=temp;cout<<temp<<"\n";}in>>temp;x=temp;cout<<"x="<<temp<<endl;in.close();}void main(){double t,y=0;read();for(int k=0;k<n;k++){t=1;for(int i=0;i<n;i++)if(k!=i)t=(x-X[i])/(X[k]-X[i])*t;y=y+t*Y[k];}cout<<"y="<<y<<endl;}1-2.埃特金逐步插值法（同上）#include<fstream.h>#include<stdlib.h>double *X,*Y; //插值节点（x，y）double x; //带估点int n; //N个插值点，插n-1次。
数值分析-源代码
一,三次样条插值#include<iostream.h>#include<stdlib.h>void ready(double*&A,double*&B,double*h,double*y,int m,int n,double m0=0,double mn=0){A=new double[n];B=new double[n];if(m==2){A[0]=0;B[0]=2*m0;A[n]=1;B[n]=2*mn;}else {A[0]=1;B[0]=3*(y[1]-y[0])/h[0];A[n]=0;B[n]=3*(y[n]-y[n-1])/h[n-1];}for(int i=1;i<n;i++){A[i]=h[i-1]/(h[i-1]+h[i]);B[i]=3*((1-A[i])*(y[i]-y[i-1])/h[i-1]+A[i]*(y[i+1]-y[i])/h[i]);}}void calculate1(double*&a,double*&b,double*A,double*B,int n){a=new double[n];b=new double[n];a[0]=-A[0]/2;b[0]=B[0]/2;for(int i=1;i<=n;i++){a[i]=-A[i]/(2+(1-A[i])*a[i-1]);b[i]=(B[i]-(1-A[i])*b[i-1])/(2+(1-A[i])*a[i-1]);}}void calculate2(double*&m,double*a,double*b,int n){m=new double[n+1];m[n+1]=0;for(int i=n;i>=0;i--)m[i]=m[i+1]*a[i]+b[i];}void main(){int e,t,w,v;out:cout<<"求三次样条插值函数:"<<endl<<"请输入函数表中x,y的对数:"<<endl;cin>>e;double*x,*y,*h,*A,*B,*a,*b,*m,s0,sn,u,s;cout<<"请输入x的一系列值:"<<endl;x=new double[e];for(int i=0;i<e;i++)cin>>x[i];cout<<"请输入y的一系列值:"<<endl;y=new double[e];for(int f=0;f<e;f++)cin>>y[f];h=new double[e];for(int d=0;d<e;d++)h[d]=x[d+1]-x[d];loop:cout<<"边界条件选择: 1.提供了s\"(x1)=0,s\"(x2)=0; 2.提供了s'(x0)=m0,s'(xn)=mn (m0,mn!=0)"<<endl;cin>>t;if(t!=1&&t!=2){cout<<"输入信息错误!请重新输入!"<<endl;goto loop;} if(t==1)ready(A,B,h,y,t,e);else {cout<<"请输入s'(x0),s'(xn)的值:"<<endl;cin>>s0>>sn;ready(A,B,h,y,t,e,s0,sn);}calculate1(a,b,A,B,e);calculate2(m,a,b,e);in:cout<<"请输入所需要计算的x的值: ("<<x[0]<<"<x<"<<x[e-1]<<")"<<endl;cin>>u;if(u<x[0]||u>x[e-1]){cout<<"输入数据错误!"<<endl;goto in;}for(int c=0;c<e;c++){if(u>x[c]&&u<x[c+1])w=c;}s=(1+2*(u-x[w])/(x[w+1]-x[w]))*((u-x[w+1])*(u-x[w+1])/((x[w]-x[w+1])*(x [w]-x[w+1])))*y[w]+(1+2*(u-x[w+1])/(x[w]-x[w+1]))*((u-x[w])*(u-x[w])/((x[w+1]-x[w])*(x[w+1 ]-x[w])))*y[w+1]+(u-x[w])*((u-x[w+1])*(u-x[w+1])/((x[w]-x[w+1])*(x[w]-x[w+1])))*m[w]+ (u-x[w+1])*((u-x[w])*(u-x[w])/((x[w+1]-x[w])*(x[w+1]-x[w])))*m[w+1]; cout<<"f("<<u<<")的近似值为: "<<endl<<s<<endl;In:cout<<"请选择:1,继续进行新的操作计算2,退出不计算"<<endl; cin>>v;if(v!=1&&v!=2){cout<<"输入错误!请重新输入!"<<endl;goto In;}if(v==1){system("cls");goto out;}}运行结果:根据课本52页的数据输入,得:第9题:f(78.3)=3.00304第10题:f(0.35)=0.591591二,自动选取步长梯形法#include<iostream.h>#include<stdlib.h>#include<math.h>double f(double x){return 2/(1+x*x);}void main(){double a,b,e,h,T0,T1,S;out:cout<<"请输入区间上下界:a和b"<<endl; cin>>a>>b;cout<<"请输入给定的精度e: "<<endl; cin>>e;h=(b-a)/2;T1=(f(a)+f(b))*h;int n=1,t,v;do{T0=T1;S=0;for(int k=1;k<n;k++)S=S+f(a+(2*k-1)*h/n);T1=T0/2+S*h/n;t=n;if(fabs(T1-T0)>=3*e)n=2*n;}while(t==n);cout<<"该积分的近似值为:"<<T1<<endl;In:cout<<"请选择:1,继续进行新的操作计算2,退出不计算"<<endl; cin>>v;if(v!=1&&v!=2){cout<<"输入错误!请重新输入!"<<endl;goto In;}if(v==1){system("cls");goto out;}}运行结果:根据课本97页第7题的数据输入,得结果为:0.75三,Romberg求积分法#include<iostream>#include<math.h>#include<stdlib.h>using namespace std;double f(double x){return sqrt(x);}void main(){double a,b,e,T[20][20],S=0,s;int k=1,v;out:cout<<"请输入区间上下界:a和b"<<endl;cin>>a>>b;cout<<"请输入给定的精度e: "<<endl;cin>>e;T[0][0]=(b-a)*(f(a)+f(b))/2;do{for(int i=1;double(i)<=pow(2,double(k))-1;i++){S=S+f(a+(2*i-1)*(b-a)/pow(2,double(k)));}T[0][k]=0.5*(T[0][k-1]+(b-a)*S/pow(2,double(k)-1));for(int m=1;m<=k;m++){T[m][k-m]=(T[m-1][k-m+1]*pow(4,double(m))-T[m-1][k-m])/(pow(4,dou ble(m))-1);}s=k;if(fabs(T[k][0]-T[k-1][0])>=e)k=k+1;}while(s!=k);cout<<"该积分的近似值为:"<<endl<<T[k][0]<<endl;In:cout<<"请选择:1,继续进行新的操作计算2,退出不计算"<<endl; cin>>v;if(v!=1&&v!=2){cout<<"输入错误!请重新输入!"<<endl;goto In;}if(v==1){system("cls");goto out;}}运行结果:根据课本97页第8题的数据输入(输入a=1,b=9,e=0.000001),得结果为:92.121四,列主元高斯消去法#include<iostream>#include<stdio.h>#include<math.h>using namespace std;void main(){int n,v;out:cout<<"请输入系数矩阵A的行列数n:"<<endl;cin>>n;double A[10][10],*b=new double[n],*x=new double[n],e,T;int i1; cout<<"请输入系数矩阵A的值: "<<endl;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++)cin>>A[i][j];}cout<<"系数矩阵A为:"<<endl;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++)cout<<A[i][j]<<" ";cout<<endl;}cout<<"请输入方程组的右端项b:"<<endl;for(int i=1;i<=n;i++)cin>>b[i];cout<<"右端项b为:"<<endl;for(int i=1;i<=n;i++)cout<<b[i]<<" ";cout<<endl;cout<<"请输入精度值e:"<<endl;cin>>e;for(int k=1;k<n;k++){double max=fabs(A[k][k]);for(int i=k;i<=n;i++){if(fabs(A[i][k])>max){max=fabs(A[i][k]);i1=i;} }T=max;if(T<e){cout<<"求解失败!"<<endl;goto stop;}if(i1!=k){ double temp1;for(int k1=1;k1<=n;k1++){double temp;temp=A[i1][k1];A[i1][k1]=A[k][k1];A[k][k1]=temp;}temp1=b[i1];b[i1]=b[k];b[k]=temp1;}for(int i=k+1;i<=n;i++){T=A[i][k]/A[k][k];b[i]=b[i]-T*b[k];for(int j=k+1;j<=n;j++)A[i][j]=A[i][j]-T*A[k][j];}}if(fabs(A[n][n])<=e){cout<<"求解失败!"<<endl;goto stop;} x[n]=b[n]/A[n][n];for(int i=n-1;i>=1;i--){double sum=0;for(int j=i+1;j<=n;j++)sum+=A[i][j]*x[j];x[i]=(b[i]-sum)/A[i][i];}for(int i=1;i<=n;i++)cout<<x[i]<<" ";cout<<endl;stop:;In:cout<<"请选择:1,继续进行新的操作计算2,退出不计算"<<endl; cin>>v;if(v!=1&&v!=2){cout<<"输入错误!请重新输入!"<<endl;goto In;}if(v==1){system("cls");goto out;}}运行结果:根据课本126页第1题的数据输入,系数矩阵为2 3 5 右端项为 5 精度为0.00001;3 4 7 61 3 3 5结果为-412五,迭代法--Jacobi迭代与Seidel迭代#include<stdio.h>#include<math.h>#include<iostream>using namespace std;void Jacobi();void Seidel();void main(){int w,v;out:loop:cout<<"请选择所要的迭代方式:1,Jacobi迭代2,Seidel迭代"<<endl;cin>>w;switch(w){case 1:Jacobi();break;case 2:Jacobi();break;default:{cout<<"输入数据错误!请重新输入!"<<endl;goto loop;} }In:cout<<"请选择:1,继续进行新的操作计算2,退出不计算"<<endl; cin>>v;if(v!=1&&v!=2){cout<<"输入错误!请重新输入!"<<endl;goto In;}if(v==1){system("cls");goto out;}}void Jacobi(){int n,k;cout<<"请输入系数矩阵A的行数n:"<<endl;cin>>n;double A[10][10],*b=new double[n],*x=new double[n],*y=new double[n],*g=new double[n],e,T;int M;cout<<"请输入系数矩阵A的值: "<<endl;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++)cin>>A[i][j];}cout<<"系数矩阵A为:"<<endl;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++)cout<<A[i][j]<<" ";cout<<endl;}cout<<"请输入方程组的右端项b:"<<endl;for(int i=1;i<=n;i++)cin>>b[i];cout<<"右端项b为:"<<endl;for(int i=1;i<=n;i++)cout<<b[i]<<" ";cout<<endl;cout<<"请输入容许误差e:"<<endl;cin>>e;cout<<"请输入容许最大迭代次数M:"<<endl;cin>>M;cout<<"请输入初始向量Y:"<<endl;for(int i=1;i<=n;i++)cin>>y[i];k=1;//--------形成迭代矩阵-------for(int i=1;i<=n;i++){if(fabs(A[i][i])<e){cout<<"求解失败!"<<endl;goto stop;}T=A[i][i];for(int j=1;j<=n;j++)A[i][j]=-A[i][j]/T;A[i][i]=0;g[i]=b[i]/T;}//--------迭代--------loop:double sum=0;for(int i=1;i<=n;i++){for(int j=1;j<=n&&j!=i;j++)sum+=(A[i][j]*y[j]);x[i]=g[i]+sum;}double fanshu=0;for(int i=1;i<=n;i++){fanshu+=fabs(x[i]-y[i]);}if(fanshu<e){for(int i=1;i<=n;i++)cout<<x[i]<<" ";cout<<endl<<k<<endl;}else if(k<M){k=k+1;for(int i=1;i<=n;i++)y[i]=x[i];goto loop;}else cout<<"求解失败!"<<endl;stop:;}void Seidel(){int n,k;cout<<"请输入系数矩阵A的行数n:"<<endl;cin>>n;double A[10][10],*b=new double[n],*x=new double[n],*y=new double[n],*g=new double[n],e,T;int M;cout<<"请输入系数矩阵A的值: "<<endl;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++)cin>>A[i][j];}cout<<"系数矩阵A为:"<<endl;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++)cout<<A[i][j]<<" ";cout<<endl;}cout<<"请输入方程组的右端项b:"<<endl;for(int i=1;i<=n;i++)cin>>b[i];cout<<"右端项b为:"<<endl;for(int i=1;i<=n;i++)cout<<b[i]<<" ";cout<<endl;cout<<"请输入容许误差e:"<<endl;cin>>e;cout<<"请输入容许最大迭代次数M:"<<endl;cin>>M;cout<<"请输入初始向量Y:"<<endl;for(int i=1;i<=n;i++)cin>>y[i];k=1;for(int i=1;i<=n;i++)x[i]=y[i];//--------形成迭代矩阵-------for(int i=1;i<=n;i++){if(fabs(A[i][i])<e){cout<<"求解失败!"<<endl;goto stop;}T=A[i][i];for(int j=1;j<=n;j++)A[i][j]=-A[i][j]/T;A[i][i]=0;g[i]=b[i]/T;}//--------迭代--------loop:double sum=0;for(int i=1;i<=n;i++){for(int j=1;j<=n&&j!=i;j++)sum+=(A[i][j]*x[j]);x[i]=g[i]+sum;}double fanshu=0;for(int i=1;i<=n;i++){fanshu+=fabs(x[i]-y[i]);}if(fanshu<e){for(int i=1;i<=n;i++)cout<<x[i]<<" ";cout<<endl<<k<<endl;}else if(k<M){k=k+1;for(int i=1;i<=n;i++)y[i]=x[i];goto loop;} else cout<<"求解失败!"<<endl;stop:;}运行结果:根据课本139页第1题的数据输入,简单迭代:x1=1.2,x2=1.4,x3=1.6,x4=0.79999Seidel迭代:x1=1.2,x2=1.4,x3=1.6,x4=0.79999。
数值分析全代码
for(k=0;k<=n-2;k++)
{
for(i=k+1;i<=n-1;i++)
{
m[i][k]=a[i][k]/a[k][k];
for(j=k+1;j<=n-1;j++)
{
for(i=0;i<n;i++)
{sum=0;
max=0;
for(j=0;j<n;j++)
{
if(i!=j)
sum=sum+a[i][j]*x[j];
}
y[i]=(b[i]-sum)/a[i][i];
c=fabs(x[i]-y[i]);
if(c>max)
max=c;
}
if(max<eps)
{printf("k=%d\n",k);
#include<math.h>
#define eps 0.0001
void main()
{
int n=3,i,j,k,n0=100;
double a[3][3]={10,-1,0,-1,10,-2,0,-2,5},b[3]={9,-5,12},t,x[3],y[3],sum=0.0,e=0.0;
for(i=0;i<n;i++)
{
sum=sum+a[i][j]*x[j];
}
x[i]=(b[i]-sum)/a[i][i];
printf("%f\n",x[i]);
数值分析数学实验上机实验编程matlab源代码
Newton法及改进的Newton法源程序：clear%%%% 输入函数f=input('请输入需要求解函数>>','s')%%%求解f(x)的导数df=diff(f);%%%改进常数或重根数miu=2;%%%初始值x0x0=input('input initial value x0>>');k=0;%迭代次数max=100;%最大迭代次数R=eval(subs(f,'x0','x'));%求解f(x0)，以确定初值x0时否就是解while (abs(R)>1e-8)x1=x0-miu*eval(subs(f,'x0','x'))/eval(subs(df,'x0','x'));R=x1-x0;x0=x1;k=k+1;if (eval(subs(f,'x0','x'))<1e-10);breakendif k>max;%如果迭代次数大于给定值，认为迭代不收敛，重新输入初值ss=input('maybe result is error,choose a new x0,y/n?>>','s');ifstrcmp(ss,'y')x0=input('input initial value x0>>');k=0; elsebreakendendendk;%给出迭代次数x=x0；%给出解Gauss消元法源程序：cleara=input('输入系数阵：>>\n')b=input('输入列阵b：>>\n')n=length(b);A=[a b]x=zeros(n,1);%%%函数主体Yipusilong-1；for k=1:n-1;%%%是否进行主元选取if abs(A(k,k))<yipusilong;%事先给定的认为有必要选主元的小数yzhuyuan=1;elseyzhuyuan=0;endifyzhuyuan;%%%%选主元t=A(k,k);for r=k+1:n;if abs(A(r,k))>abs(t)p=r;else p=k;endend%%%交换元素if p~=k; for q=k:n+1;s=A(k,q);A(k,q)=A(p,q);A(p,q)=s;endendend%%%判断系数矩阵是否奇异或病态非常严重if abs(A(k,k))<yipusilongdisp(‘矩阵奇异，解可能不正确’)end%%%%计算消元,得三角阵for r=k+1:n;m=A(r,k)/A(k,k);for q=k:n+1;A(r,q)=A(r,q)-A(k,q)*m;endendend%%%%求解xx(n)=A(n,n+1)/A(n,n);for k=n-1:-1:1;s=0;for r=k+1:n;s=s+A(k,r)*x(r);endt=(A(k,n+1)-s) x(k)=(A(k,n+1)-s)/A(k,k) end Lagrange插值源程序：n=input('将区间分为的等份数输入：\n');s=[-1+2/n*[0:n]];%%%给定的定点，Rf为给定的函数x=-1:0.01:1;f=0;for q=1:n+1;l=1;%求插值基函数for k=1:n+1;if k~=q;l=l.*(x-s(k))./(s(q)-s(k));elsel=l;endendf=f+Rf(s(q))*l;%求插值函数endplot(x,f,'r')%作出插值函数曲线grid onhold on分段线性插值源程序clearn=input('将区间分为的等份数输入：\n');s=[-1+2/n*[0:n]];%%%给定的定点，Rf为给定的函数m=0;hh=0.001;for x=-1:hh:1;ff=0;for k=1:n+1;%%%求插值基函数switch kcase 1if x<=s(2);l=(x-s(2))./(s(1)-s(2));elsel=0;endcase n+1if x>s(n);l=(x-s(n))./(s(n+1)-s(n)); elsel=0;endotherwiseif x>=s(k-1)&x<=s(k);l=(x-s(k-1))./(s(k)-s(k-1)); else if x>=s(k)&x<=s(k+1);l=(x-s(k+1))./(s(k)-s(k+1)); elsel=0;endendendff=ff+Rf(s(k))*l;%%求插值函数值endm=m+1;f(m)=ff;end%%%作出曲线x=-1:hh:1;plot(x,f,'r');grid onhold on三次样条插值源程序：（采用第一边界条件）clearn=input('将区间分为的等份数输入：\n');%%%插值区间a=-1;b=1;hh=0.001;%画图的步长s=[a+(b-a)/n*[0:n]];%%%给定的定点，Rf为给定的函数%%%%第一边界条件Rf"(-1),Rf"(1)v=5000*1/(1+25*a*a)^3-50/(1+25*a*a)^4;for k=1:n;%取出节点间距h(k)=s(k+1)-s(k);endfor k=1:n-1;%求出系数向量lamuda,miula(k)=h(k+1)/(h(k+1)+h(k));miu(k)=1-la(k);end%%%%赋值系数矩阵Afor k=1:n-1;for p=1:n-1;switch pcase kA(k,p)=2;case k-1A(k,p)=miu(p+1);case k+1A(k,p)=la(p-1);otherwiseA(k,p)=0;endendend%%%%求出d阵for k=1:n-1;switch kcase 1d(k)=6*f2c([s(k) s(k+1) s(k+2)])-miu(k)*v; case n-1d(k)=6*f2c([s(k) s(k+1) s(k+2)])-la(k)*v; otherwised(k)=6*f2c([s(k) s(k+1) s(k+2)]);endend %%%%求解M阵M=A\d';M=[v;M;v];%%%%m=0;f=0;for x=a:hh:b;if x==a;p=1;elsep=ceil((x-s(1))/((b-a)/n));endff1=0;ff2=0;ff3=0;ff4=0;m=m+1;ff1=1/h(p)*(s(p+1)-x)^3*M(p)/6;ff2=1/h(p)*(x-s(p))^3*M(p+1)/6;ff3=((Rf(s(p+1))-Rf(s(p)))/h(p)-h(p)*(M(p+1)-M(p))/6)*(x-s(p));ff4=Rf(s(p))-M(p)*h(p)*h(p)/6;f(m)=ff1+ff2+ff3+ff4 ;end%%%作出插值图形x=a:hh:b;plot(x,f,'k')hold on grid on 多项式最小二乘法源程序clear%%%给定测量数据点(s,f)s=[3 4 5 6 7 8 9];f=[2.01 2.98 3.50 5.02 5.47 6.02 7.05];%%%计算给定的数据点的数目n=length(f);%%%给定需要拟合的数据的最高次多项式的次数m=10;%%%程序主体for k=0:m;g=zeros(1,m+1);for j=0:m;t=0;for i=1:n;%计算内积(fai(si),fai(si))t=t+fai(s(i),j)*fai(s(i),k);endg(j+1)=t;endA(k+1,:)=g;%法方程的系数矩阵t=0;for i=1:n;%计算内积(f(si),fai(si))t=t+f(i)*fai(s(i),k);endb(k+1,1)=t;enda=A\b%求出多项式系数x=[s(1):0.01:s(n)]';y=0;fori=0:m;y=y+a(i+1)*fai(x,i);endplot(x,y)%作出拟合成的多项式的曲线grid onhold onplot(s,f,'rx') %在上图中标记给定的点表中，L10(x)为Lagrang e插值的10次多项式，S10(x)，S40(x)分别代表n=10，40的三次样条插值函数，X10(x)，X40(x)分别代表n=10，40的线性分段插值函数。
一些经典的数值分析(matlab程序)
1、牛顿迭代法此方法一般用来求函数的根，速度比较快，程序比较简单。
文件newton.m 内容如下%n表示迭代次数，ret为返回值function ret=newton(n)format longy=inline('x^3+10*x-20');z=inline('3*x^2+10');x0=1.5;for i=1:na=y(x0);b=z(x0);x1=x0-a/b;x0=x1;endret=x12、复合Simpson求积分很简单的一种求定积分的方法，将求积区间分成很多小的曲边梯形，再累加。
文件mulsimpson.m内容如下%复化simpson求积分%a,b表示区间，n表示区间数function ret=mulsimpson(a,b,n)h=(b-a)/n;detsum=0;for i=1:n-1xk=a+i*h;detsum=detsum+fun(xk);endret=h*(fun(a)+fun(b)+2*detsum)/2;%内建函数function z=fun(x)z=exp(-x*x);3、4阶龙格-库塔法求积分此方法速度较快，编程较方便文件step4Runge.m内容如下%4阶Runge-Kutta 法%a,b为积分区间，N为划分数目，y0为初值，函数由fun定义function ret=step4Runge(a,b,N,y0)format longh=(b-a)/N;n=1;x0=a;for n=1:Nx=x0+h;k1=fun(x0,y0);k2=fun(x0+h/2,y0+h*k1/2);k3=fun(x0+h/2,y0+h*k2/2);k4=fun(x0+h,y0+h*k3);y=y0+h*(k1+2*(k2+k3)+k4)/6;x0=xy0=yend%积分函数定义function z=fun(x,y)z=1+y*y;4、Gauss_Seidel迭代解线性方程相比消元法，编程较为容易文件Gauss_Seidel.m内容如下%此函数演示高斯-赛德尔迭代%a表示系数矩阵，b表示值矩阵，n表示系数矩阵阶数，M表示迭代次数%注意b为列向量function y=Gauss_Seidel(a,b,n,M)format longx0=[0;0;0];for k=1:Mfor i=1:ns=0;t=x0(i);for j=1:nif j~=is=s+a(i,j)*x0(j);endendx0(i)=(b(i)-s)/a(i,i);endendy=x0;5、高斯列主消元法此方法为解线性方程组常用方法，先化简增广矩阵，然后回代求解。
数值分析python代码
from math import logfrom math import sin,cosimport random#Q1def bifection(f,a,b,error,n):"""f stands for the function we want to approximate to, while a and b are the left and right number of the interval where the f effects in, and n means for thetime we have tried to make a approximation""""""bifection(f,0,1,10**-6,20)"""FA = f(a)for i in range(n):p = a + (b-a)/2FP = f(p)if FP == 0 or (b - a)/2 < error:return pif FA * FP > 0:a = pelse:b = pdef iteration(f,p0,n,error):"""iteration(f,0.5,22,10**(-6))"""p = 0for i in range(n):p = (2.71828)**(-p0)if abs(p - p0) < error:return pelse:p0 = pdef fastitertion(f,p0,n,error):"""fastitertion(f,0.5,20,10**-6)"""p = 0q0 = p0for i in range(n):p = (2.71828)**(-p0)q = p + g(p)/(1 - g(p))*(p - q0)if abs(q - q0) < error:return qelse:p0 = pq0 = qdef newton(f,g,p0,n,error):"""newton(f,0.5,3,10**-6)"""for i in range(n):p = p0 - f(p0)/g(p0)if abs(p - p0) < error:return pelse:p0 = pdef secant(f,p0,p1,n,error):"""secont(f,0.5,1,6,10**-6)"""for i in range(n):p = p0 - f(p0)*(p1 - p0)/(f(p1) - f(p0)) if abs(p - p0) < error:return pelse:p0 = pdef fastsecant(f,p0,p1,n,error):"""secont(f,0.5,1,4,10**-6)"""for i in range(n):p = p0 - f(p0)*(p1 - p0)/(f(p1) - f(p0)) if abs(p - p0) < error:return pelse:p1 = p0p0 = pdef f(x):return x - pow(2.71828, - x)def g(x):return 1 + pow(2.71828, - x)#Q2def q2_eqp(x):return 2*(2.71828)**(-x) - sin(x)def q2_eqq(x):return -2*(2.71828)**(-x) - cos(x)def q2():return newton(q2_eqp,q2_eqq,3,20,10**-12) def q3():return newton(q2_eqp,q2_eqq,0,20,10**-12)#Q3def F(x,y):return x**2 + y**2 - 1def G(x,y):return x**3 - ydef f1(x):return 2 * xdef f2(x):return 2 * xdef g1(x):return 3 * x**2def g2(x):return - 1def newtonforsystem(p,q,x0,y0,error):"""newtonforsystem(F,G,0.8,0.6,error)[[1.0, 0.0, 0.83013698630137], [0.0, 1.0, 0.5698630136986301]]"""lst1 = [[f1(x0),f2(x0)],[g1(x0),g2(x0)]]lst2 = [- p(x0,y0) + f1(x0) * x0 + f2(x0) * y0 , - q(x0,y0) + g1(x0) * x0 + g2(x0) * y0]lst = [[f1(x0),f2(x0),lst2[0]],[g1(x0),g2(x0),lst2[1]]]return gaussjordaneli(lst)#Q4def itermeth(lst1,lst2,n,x_lst):"""this is a programme which stands for the simple iteration method to solve the system of linear equation. where the lst is analogy for the matrix.lst1 = [[8,1,2],[4,-10,-1],[-3,-2,5]],lst2 = [12,21,16],error = 10**-5answer: [1.0, -2.0, 3.0]"""temp_lst = [0 for x in range(len(lst1[0]))]temp = x_lst[0]for i in range(len(lst1)):lst2[i] = lst2[i] / lst1[i][i]lst1[i] = [- item / lst1[i][i] for item in lst1[i]]for i in range(len(lst1)):lst1[i][i] = 0lst1[i].append(lst2[i])for n in range(n):for i in range(len(x_lst)):temp_lst[i] = sum([x_lst[j] * lst1[i][j] for j in range(len(lst1[i]) - 1)]) + lst1[i][-1]x_lst = temp_lstreturn temp_lstdef seidelitermeth(lst1,lst2,n,x_lst):x_lst = itermeth(lst1,lst2,1)temp_lst = x_lst[0] + x_lst[1:]x_lst2 = itermeth(lst1,lst2,n,temp_lst)temp_lst = x_lst[0] + x_lst2[1] + x_lst[-1]return itermeth(lst1,lst2,n,temp_lst)#Q5#target matrix is [[2,3,-1,4,1,-5],[3,0,4,-1,-5,0],[5,-1,1,-4,1,8],[-1,3,4,5,-3,4],[4,-4,3,7,5,1]]def gaussjordaneli(group):""">>> group = [[2,3,-1,4,1,-5],[3,0,4,-1,-5,0],[5,-1,1,-4,1,8],[-1,3,4,5,-3,4],[4,-4,3,7,5,1]] >>> gausseli(group)[1 0 0 0 0 -1][0 1 0 0 0 2][0 0 1 0 0 4][0 0 0 1 0 -2][0 0 0 0 1 3]"""for i in range(len(group)):first = group[i][i]for j in range(i,len(group[0])):group[i][j] = group[i][j] / firstfor temp in range(len(group)):if temp != i:second = group[temp][i]for j in range(i,len(group[0])):group[temp][j] -= group[i][j] * secondreturn groupdef gausseli(group):"""[[1.0, 1.5, -0.5, 2.0, 0.5, -2.5][-0.0, 1.0, -1.2, 1.5, 1.4, -1.6][-0.0, -0.0, 1.0, 0.1, -1.5, -0.9][-0.0, -0.0, -0.0, 1.0, -5.0, -17][ 0.0, 0.0, 0.0, 0.0, 1.0, 3.0]"""for i in range(len(group)):first = group[i][i]for j in range(len(group[0])):group[i][j] = group[i][j] / firstfor temp in range(i + 1,len(group)):second = group[temp][i]for j in range(i,len(group[0])):group[temp][j] -= group[i][j] * secondreturn groupdef jordan(group):"""[[2, 3, -1, 4, 1, -5][0.0, -4.5, 5.5, -7.0, -6.5, 7.5][0.0, 0.0, -6.9, -0.8, 10.8, 6.3][0.0, 0.0, 0.0, -1.0, 5.1, 17.3][0.0, 0.0, 0.0, 0.0, 83.0 249.0]"""for i in range(len(group)):first = group[i][i]for temp in range(i + 1,len(group)):second = group[temp][i]for j in range(i,len(group[0])):group[temp][j] -= group[i][j] / first * secondreturn group#fittingdef fitting(lst,f,g):"""f and g is the function which translate the original function into linear"""i = 0total = 0lst[0] = [g(item) for item in lst[0]]lst[1] = [f(item) for item in lst[1]]x_bar = sum(lst[0])/len(lst[0])y_bar = sum(lst[1])/len(lst[1])while i < len(lst[0]):total = total + lst[0][i] * lst[1][i]i = i + 1b = (total - len(lst[0]) * x_bar * y_bar)/(sum([x**2 for x in lst[0]]) - len(lst[0]) * (x_bar**2))a = y_bar -b * x_barreturn (b,a)#Q6#ln(y) = bx + ln(a)""">>>lst1 = [[1,2,3,4,5,6,7,8],[15.3,20.5,27.4,36.6,49.1,65.6,87.8,117.6]] >>>a = fitting(lst1, log, lambda x: x)a = (0.29121601623818705, 2.436859706328185)"""#a = e^(a)#y = 11.4375e^(0.291x)#Q7#1/y = b/t + a""">>>lst2 = [[1,2,3,4,5,6,7,8],[4.00,6.40,8.00,8.80,9.22,9.50,9.70,9.86]] >>>b = fitting(lst2,lambda y: 1/y, lambda x: 1/x)b = (0.17160826360564463, 0.07458941352083057)"""#y = t/(0.075t + 0.171)#Q8#y = a + bx^3"""vender = [-5,-4,-3,-2,-1,0,1,2,3,4,5]beta = [-12.9,-5.95,-1.76,0.42,1.20,1.34,1.43,2.25,4.38,8.61,15.6]lst = [vender,beta]b = fitting(lst,lambda y: y, lambda x: x**3)(0.11, 1.33)"""#y = 1.33 + 0.11x^3def matrix_trans(matrix):"""[[a,b] [[a,c][c,d]] ----> [b,d]]"""return [[row[i] for row in matrix] for i in range(len(matrix[0]))]def matrix_inver(matrix):for temp in range(len(matrix)):matrix[temp].extend([0 for x in range(len(matrix))])for i in range(len(matrix)):matrix[i][len(matrix) + i] = 1inverse = gaussjordaneli(matrix)for i in range(len(inverse)):inverse[i] = inverse[i][len(inverse):]return inversedef matrix_mul(A, B):result = [[0] * len(B[0]) for i in range(len(A))]for i in range(len(A)):for j in range(len(B[0])):for k in range(len(B)):result[i][j] += A[i][k] * B[k][j]return resultdef polyfitting(vender,beta):"""vender = [[1,0,125],[1,0,64],[1,0,27],[1,0,8],[1,0,1],[1,0,0],[1,0,1],[1,0,8],[1,0,27],[1 ,0,64],[1,0,125]]beta = [-12.9,-5.95,-1.76,0.42,1.20,1.34,1.43,2.25,4.38,8.61,15.6]"""new_beta = [[beta[i]] for i in range(len(beta))]temp = matrix_mul(matrix_inver(matrix_mul(matrix_trans(vender),vender)),matrix_trans(ve nder))result, i = [] , 0while i < len(temp):result.append(sum([temp[i][j] * beta[j] for j in range(len(beta))]))i += 1return result#Q9 need to finish the gausseli first"""p = [[1,1,1],[1,2,4],[1,3,9],[1,4,16]]q = [4,10,18,26][-1.5, 4.9, 0.5]"""#Q10def lagrangeitp(lst,x):""">>>lst =[[1.0,1.5,2.0,3.5,5.5,6.0],[1.0,3.375,8.0,42.875,166.375,216.0]]>>>lagrangeitp(lst,2.5)15.625>>>lagrangeitp(lst,4.5)91.125"""return sum([(la_multool(lst[0],i,x) * lst[1][i]) for i in range(len(lst[0]))])def la_multool(lst,i,x):k = 0total = 1while k < len(lst):if k != i:total = total * (x - lst[k])/(lst[i] - lst[k])k = k + 1return total#newton interpolationdef newtonitp(lst,x):return sum([newton_multool(lst,0,i) * multool(lst[0],i - 1, x) for i in range(1,len(lst[0]))]) + lst[1][0]def multool(lst,i,x):if i == 0:return x - lst[0]else:return (x - lst[i]) * multool(lst, i-1, x)def newton_multool(lst,i,j):"""lst is a list which has rank 2,the first is x whilethe second holds for y."""if i > j:return 0if i == j :return lst[1][i]else:return (newton_multool(lst,i+1,j) - newton_multool(lst,i,j - 1))/(lst[0][j] - lst[0][i])#Q11""">>>lst = [[0.4,0.5,0.6,0.7,0.8,0.9],[-0.916291,-0.693147,-0.510826,-0.357765,-0.223144,-0 .105361]]>>>newtonitp(lst,0.78)-0.24887097759999988>>>log(0.78)-0.2484613592984996error = 0.0016"""#Q12""">>>lst = [[[0.785, 0.873, 0.959, 1.047]],[0.7071,0.7660,0.8192,0.8660]]>>>newtonitp(lst,0.907)0.7071>>>sin(0.907)0.7876error = 0.1"""#Q13""">>>lst = [[0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9],[0.4794,0.6442,0.7833,0.8912,0.9636,0.9975,0. 9917,0.9463]]>>>q1 = lagrangeitp([lst[0][0:3],lst[1][0:3]],0.47)0.452463375>>>sin(0.47)0.452886285error = 0.0008>>>q2 = lagrangeitp([lst[0][4:7],lst[1][4:7]],1.6)0.9995625>>>sin(1.6)0.9995736error = 10^(-5)"""class integral(object):def __init__(self,function,a,b):self.f = functionself.a = aself.b = bdef simpson(self,n):h = (self.b - self.a)return h * (self.f(self.a) + 4 * self.f(1/2 *(self.a + self.b)) + self.f(self.b)) / 6def fixed_simpson(self,n):h = (self.b - self.a)/ nresult = self.f(self.a) + self.f(self.b)for k in range(n):result += 4 * self.f(self.a + (k - 1/2)*h) + 2 * self.f(self.a + k * h) return h * result / 6def fixed_cotes(self,n):h = (self.b - self.a)/ nresult = (self.f(self.a) + self.f(self.b)) * 7for k in range(n):result += 32 * self.f(self.a + (k - 1/4)* h) + 12 * self.f(self.a + (k - 1/2) * h) +\32 * self.f(self.a + (k - 3/4)* h) + 14 * self.f(self.a + k * h)return h * result / 90def romberg(self,n):R = [[0 for x in range(n)] for x in range(n)]h = (self.b - self.a)/ nR[0][0] = (self.f(self.a) + self.f(self.b)) * h / 2for i in range(n - 1):for j in range(i):R[i][j] = self.T(i,j,n)return Rdef T(self,i,k,n):if i == 0 and k == 0:return (self.f(self.a) + self.f(self.b)) * (self.b - self.a)/ (n * 2) if i == 0 and k != 0:return 1/2 * self.T(0,k - 1,n) + (self.b - self.a) / (2 ** k) * sum([self.f(self.a + (2*i + 1)*(self.b - self.a)/(2 ** k)) for i in range(2 ** (k - 1) - 1)])else:return (4 ** i * self.T(i - 1,k+1,n) - self.T(i - 1,k,n))/(4**i - 1)#Q14""">>> E1 = integral(lambda x:(2.71828)**(-x),0,1)>>> E1.simpson(90)0.6323338572410494>>> original = 0.6321203113733684"""#Q15""">>>s1 = integral(lambda x: sin(x)/x ,0 + 1/64,5)>>> s1.fixed_simpson(64)1.622325**********>>>s1 = integral(lambda x: sin(x)/x ,0 + 1/128,5)>>> s1.fixed_simpson(128)1.5861535544830339>>>s1 = integral(lambda x: sin(x)/x ,0 + 1/256,5)>>> s1.fixed_simpson(256)1.5680481231793797"""#Q16""">>>y1 = integral(lambda x : 1/x , 1, 3)>>>y1.romberg(4)"""class diff_equation(object):def __init__(self,function,initial):"""function = lambda x : lambda y :~"""self.f = functionself.initial = initialdef Euler(self,a,b,h):N = (b - a)/htemp = aresult = self.initialfor i in range(int(N)):print(result)result += h * self.f(temp + i * h, result)temp = temp + hreturn resultdef improved_euler(self,a,b,h):N = (b - a)/htemp = aresult = self.initialfor i in range(int(N)):print(result)temp1 = resultresult += h / 2 * (self.f(temp + i * h, temp1 + h * self.f(temp, temp1)) + self.f(temp,temp1))temp = temp + hreturn resultdef Runge_Kutta4(self,a,b,h):N = (b - a)/htemp = aresult = self.initialK1 = self.f(a,result)K2 = self.f(a + 1/2 * h, result + 1 / 2 * K1 * h)K3 = self.f(a + 1/2 * h, result + 1 / 2 * K2 * h)K4 = self.f(a + h, result + K3 * h)for i in range(int(N)):print(result)result += 1/6 * h *(K1 + K2 + K3 + K4)temp = temp + hK1 = self.f(temp,result)K2 = self.f(temp + 1/2 * h, result + 1 / 2 * K1 * h) K3 = self.f(temp + 1/2 * h, result + 1 / 2 * K2 * h) K4 = self.f(temp + h, result + K3 * h)return result#Q17"""0.00.20.520.8080.96161.01.01.01.01.0[0.0,0.39346913629508384,0.6321203113733684,0.7768696147177936,0.8646645346959726,0.9179148633392415,0.950212831163814,0.9698025454843657,0.9816843118309424,0.9888909698354715]"""#Q18""">>> f = lambda x,y : x + y>>> k = diff_equation(f,1)>>> k.Runge_Kutta4(0,1,0.2)11.16213333333333321.34821916444444431.59134274233837041.92947711013982562.4062051952111494>>> k.improved_euler(0.1,0.2)11.221.55242.0219282.658752163.4996776352[1.0, 1.2428051876884079, 1.5836485924986987, 2.0442361299975187,2.651079461759733]"""#Q19def diff_sys(b,interval,f0,g0):f = lambda x,y : x * (1 - x ** 2 - y ** 2) + y * (x**2 - y**2 - b)g = lambda x,y : y * (1 - x ** 2 - y ** 2) - x * (x**2 - y**2 - b)N = (interval[1] - interval[0])/0.05temp = interval[0]result = [f0,g0]for i in range(int(N)):print(result)result[0] += 0.05 * f(result[0],result[1])result[1] += 0.05 * g(result[0],result[1])temp = temp + 0.05return result#Q20def partial(h,k):"""pU/pT = [u(n+1,k) - u(n,k)]/dtp^2U/px^2 = [u(n,k+1) - 2 * u(n,k) +u(n,k - 1)]/dx^2u|t = 0 \ u(0,i) = 4(i*dx)*(1 - i*dx)u|x = 0 \ = u| x = 1\ = 0"""#Q22""">>> monte_carlo(100,lambda x:2.71828 ** x) 1.6037852>>> monte_carlo(10000,lambda x:2.71828 ** x) 1.7225740360000001"""def monte_carlo(number,f,a,b):count = 0for i in range(1,number + 1):x = random.uniform(0,a)y = random.uniform(0,b)if y < f(x):count += 1return count / number * (a * b)。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值积分
h=0.01;
x=0:h:1;
y=4./(1+x.^2);
format long
t=length(x); %数组长度
z1=sum(y(1:(t-1)))*h %矩形公式2
z2=sum(y(2:t))*h %矩形公式2
z3=trapz(x,y) %梯形公式
z4=quad('4./(1+x.^2)',0,1) %辛普森公式
方程组的解法1-雅可比迭代
X0=[0,0,0]';
A=[10,-1,-2;-1,10,-2;-1,-1,5];
b=[7.2,8.3,4.2]';
X1=A\b;
t=[];
n=10;
for k=1:n
for j=1:3
X(j)=(b(j)-A(j,[1:j-1,j+1:3])*X0([1:j-1,j+1:3]))/A(j,j); end
t=[t;X'];
X0=X;
end
disp('方程组的解')
X1
disp('迭代10步的解')
X
disp('方程组每次迭代的解')
t
方程组的解法2—高斯塞德尔迭代
X=[0,0,0]';
A=[10,-1,-2;-1,10,-2;-1,-1,5];
b=[7.2,8.3,4.2]';
X1=A\b;
t=[];
n=10;
for k=1:n
for j=1:3
X(j)=(b(j)-A(j,[1:j-1,j+1:3])*X([1:j-1,j+1:3]))/A(j,j);
end
t=[t,X];
end
disp('方程组的解')
X1
disp('迭代10步的解')
X
disp('方程组每次迭代的解')
t'
常微分方程
x0=0;
xf=1;
y0=pi/2;
[x,y]=ode23('funst',[x0,xf],y0);
yy=(x+pi/2)./cos(x);%精确解
[x,y,yy]
不要忘记funst.m也是个文件名，调用的，和这段程序放在同一文件夹下面
插值方法
x=0:1:6;
y=cos(x);
xi=0:0.25:6;
yi1=interp1(x,y,xi,'linear');
yi2=interp1(x,y,xi,'cubic');
plot(x,y,'*',xi,yi1,xi,yi2)
在MATLAB中，一维多项式插值的方法通过命令interp1(x,y,xi,method)实现，其具体的调用格式如下：
插值的方法method参数的取值和对应的含义如下：
∙linear：线性插值（linear interpolation）
∙pchip：分段三次厄米特多项式差值（piecewise cubic Hermite interpolation）。
∙cubic：三次多项式插值，与分段三次厄米特插值方法相同。
各种插值方法的比较
方法说明
linear 执行速度较快，有足够的精度，最为常用，而且为默认设置。
cubic 较慢，精度高，平滑度好，当希望得到平滑的曲线时可以使用该选项。