数值分析实验代码

合集下载

数值分析实验代码

埃特金#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;}。

数值分析实验报告程序

一、Newton插值法#include<stdio.h>#define MAX_N 20typedef struct tagPOINT{double x;double y;}POINT;int main(){int n,i,j;POINT points[MAX_N+1];double diff[MAX_N+1]; double x,tmp,newton=0;printf("\nInput n value:");scanf("%d",&n);printf("Now input the (x_i,y_i),....,%d:\n",n);for(i=0;i<=n;i++)scanf("%lf%lf",&points[i].x,&points[i].y);printf("Now input the x value:");scanf("%lf",&x);for(i=0;i<=n;i++) diff[i]=points[i].y;for(i=0;i<=n;i++){for(j=n;j>i;j--){diff[j]=(diff[j]-diff[j-1])/(points[j].x-points[j-1-i].x); }}tmp=1;newton=diff[0];for(i=0;i<n;i++){tmp=tmp*(x-points[i].x);newton=newton+tmp*diff[i+1];}printf("newton(%f)=%f\n",x,newton);return 0;}二、Lagrange插值法double Lagrange(int n,double u,double x[],double y[]) {int i,j;double l,Ln=0;for(i=0;i<=n;i++){l=1.0;for(j=0;j<=n;j++)if (i!=j)l=l*(u-x[j])/(x(i)-x[j]);Ln=Ln+l*y[i];}return(Ln);}#include<stido.h>#define N 20V oid main(void){int i,n;double x[N],y[N],u.fu;double Lagrange(int n,double u,double x[],double y[]); printf(“请输入插值多项式的次数N及插入点U：\n”); scanf(“%d,%lf”,&n,&u);printf(“请输入n+1个插值节点x,y:\n”);for(i=0;i<=n;i++)scanf(“%lf,%lf”,&x[i],&y[i]);fu=Lagrange(n,u,x,y);printf(“插入点u的近似值为%lf\n”,fu);}三、用梯形公式#include<stdio.h>#include<math.h>double f(double x){return(sqrt(pow(x,4)+1));}double T(double a,double b){double u,h;h=(b-a)/2;u=f(a)+f(b);return (h*u);}void main(){ double a,b;printf("请输入积分区域：a,b\n"); scanf("%lf,%lf",&a,&b);double T(double ,double );printf("%lf,%lf",T(a,b));四、Simpson公式求积分#include<stdio.h>#include<math.h>double f(double x){return(sqrt(2-pow(sin(x),2))); }double S(double a,double b){double u,h,c;h=(b-a)/6;c=(a+b)/2;u=f(a)+f(b)+4*f(c);return (h*u);}void main(){ double a,b;printf("请输入积分区域：a,b\n");scanf("%lf,%lf",&a,&b);double S(double ,double );printf("%lf ",S(a,b))五、复化数值积分公式的自动控制误差算法复化Simpson自动控制误差#include<math.h>double AS(double (*f)(double x),double a,double b,double e) {int i;int n=2;double h,Jn,J2n,Sn,S2n;h=(b-a)/n;Jn=(*f)((a+b)/2);Sn=h/3*((*f)(a)+(*f)(b)+4*Jn);while(1){J2n=0;for(i=1;i<=n;i++)J2n=J2n+(*f)(a+(2*i-1)*h/2);S2n=Sn/2+h/6*(4*J2n-2*Jn);if(fabs(S2n-Sn)<15*e){return S2n;break;}else{ h=h/2;n=2*n;J2n=Jn;Sn=S2n;}}}复化梯形公式double comTx(double (*f)(double),double a,double b,int n) {int i; double h,Tn;h=(b-a)/n;Tn=(*f)(a)+(*f)(b);for(i=1;i<=n-1;i++)Tn=Tn+2*(*f)(a+i*h);Tn=(h/2)*Tn;return(Tn);}复化梯形自动控制误差double A T(double (*f)(double x),double a,double b,double e) {int i; int n=10;double h,Js,Tn,T2n;h=(b-a)/n;Tn=comTx(f,a,b,n);while(1){Js=0;for(i=1;i<=n;i++)Js=Js+(*f)(a+(2*i-1)*h/2);T2n=Tn/2+h/2*Js;if(fabs(T2n-Tn)<3*e)break;else{ h=h/2;n=2*n;Tn=T2n;}}return T2n;}六、Romberg公式计算积分#include "stdio.h"#include "math.h"#define MAX 20double f(double x)//要求的积分公式{ return(lnx);}double computeTn(double (*f)(double x),double a,double b, int n)//复化梯形公式{int i;double sum=0;double h=(b-a)/n;for(i=1;i<n;i++)sum+=f(a+i*h);sum+=(f(a)+f(b))/2;return (h*sum);}double Romberg(double (*f)(double x),double a,double b,double ep)//Romberg求积分{ int j,n=1,k;double R[3][MAX];R[2][1]=computeTn(f,a,b,n);printf("%lf\n",R[2][1]);for(k=2;k<MAX;k++){ for(j=1;j<k;j++)R[1][j]=R[2][j];n*=2;R[2][1]=computeTn(f,a,b,n);printf("%lf\t",R[2][1]);for(j=2;j<=k;j++){ R[2][j]=R[2][j-1]+(R[2][j-1]-R[1][j-1])/(pow(4,j-1)-1);printf("%lf\t",R[2][j]);}printf("\n");if(fabs(R[1][k-1]-R[2][k])<ep) return(R[2][k]); }k--;if(fabs(R[1][k-1]-R[2][k])>ep){printf("Return no solved...\n");return 0;}}七、用牛顿迭代法、弦截法求非线性方程的根///用牛顿迭代法的C源程序///#include<stdio.h>#include<math.h>#define m 20double f(double x){return((x*(x*(x-7.7)+19.2)-15.3));}double f1(double x){return(x*(3*x-15.4)+19.2);}main(){int k;double x,x0,e;printf("请输入x0和允许误差e的值:\n");scanf("%lf,%lf",&x0,&e);for(k=1;k<=m;k++){x=x0-f(x0)/f1(x0);if(fabs(x-x0)<e){printf("方程的近似根为%lf",x);return NULL;} x0=x;}printf("给定的迭代次数太小.停机");}///用弦截法的C源程序///#include<stdio.h>#include<math.h>#define m 20double f(double x){return((x*(x*(x-7.7)+19.2)-15.3));}main(){int k;double x,x0,x1,f0,f1,e;printf("请输入x0,x1和允许误差e的值:\n");scanf("%lf,%lf,%lf",&x0,&x1,&e);f0=f(x0);for(k=1;k<=m;k++){f1=f(x1);x=x1-(f1*(x1-x0)/(f1-f0));if(fabs(x-x1)<e){printf("方程的近似根为%lf",x);return NULL;} f0=f1;x0=x1;x1=x;}printf("给定的迭代次数太小.停机");}八、Gauss-Seidel迭代法解线性方程组的算法#include <math.h>#include <stdio.h>#define eps 0.000001#define max 100#define N 20double norm_inf(double x[],int n){ double norm;int i;norm=fabs(x[0]);for(i=1;i<n;i++){if(fabs(x[i])>norm)norm=fabs(x[i]);}return norm;}void seidel(double a[N][N],double g[N],int n){double b[N][N]={0},x0[N],x1[N],x1_x0[N],norm,temp; int i,j,k;for(i=0;i<n;i++){g[i]=g[i]/a[i][i];for(j=0;j<n;j++){ if(i==j)continue;b[i][j]=-a[i][j]/a[i][i];}}for(i=0;i<n;i++){x0[i]=0;x1[i]=1;x1_x0[i]=x1[i]-x0[i];}k=0;norm=norm_inf(x1_x0,n); while((norm>=eps)&&(k<max)) {for(i=0;i<n;i++)x0[i]=x1[i];for(i=0;i<n;i++){temp=0;for(j=0;j<=i-1;j++)temp=temp+b[i][j]*x1[j];for(j=i+1;j<n;j++)temp=temp+b[i][j]*x0[j];x1[i]=temp+g[i];x1_x0[i]=x1[i]-x0[i];}norm=norm_inf(x1_x0,n);k++;}for(i=0;i<n;i++)printf("x[%d]=%lf\n",i,x1[i]); printf("%d times iteration.\n",k); }int main(){double b[N][N]={{10,-2,-1},{-2,10,-1},{-1,-2,5}},f[N]={0,-21,-20}; printf("\n");seidel(b,f,3);return 0;}九、用列主元消元法#include<stdio.h>#include<math.h>#define n 3void print(double a[n][n+1]);void gauss(double a[n][n+1],double x[n]){ int i,j,k;double temp,s,l;for(i=0;i<n-1;i++){ k=i;for(j=i+1;j<n;j++){ if(fabs(a[j][i])>fabs(a[k][i]))k=j;}if(k!=i)for(j=i;j<=n;j++) {temp=a[i][j];a[i][j]=a[k][j];a[k][j]=temp;}for(j=i+1;j<n;j++){ l=1.0*a[j][i]/a[i][i];for(k=0;k<n+1;k++)a[j][k]=a[j][k]-a[i][k]*l; }print(a);printf("lf");}print(a);x[n-1]=a[n-1][n]/a[n-1][n-1]; for(i=n-2;i>=0;i--){ s=0.0;for(j=i;j<n;j++){ if(j==i)continue;s+=a[i][j]*x[j];}x[i]=(a[i][n]-s)/a[i][i];}}void print(double a[n][n+1]){int i,j;for(i=0;i<n;i++){ for(j=0;j<n+1;j++)printf("%lf",a[i][j]);printf("lf");}}十、Doolittle分解法#include "iostream.h"#include "stdio.h"#define max 20void main(){int i,k,n,j,p;double x,v,y;double a[max][max],r[max][max],l[max][max]; double b[max],Y[max],X[max];printf("Input n:"); scanf("%d",&n);printf("Input a[%d][%d]:\n",n,n);for(i=1;i<n+1;i++)for(j=1;j<n+1;j++){ scanf("%lf",&x); a[i][j]=x; }printf("b[i]:");for(i=1;i<n+1;i++){ scanf("%lf",&y); b[i]=y; }for(k=1;k<n+1;k++){ for(i=k;i<n+1;i++){ v=0;for(p=1;p<k;p++)v=l[k][p]*r[p][i]+v;r[k][i]=a[k][i]-v; printf("r[%d][%d]=%f\n",k,i,r[k][i]); } for(i=k+1;i<n+1;i++){ v=0;for(p=1;p<k;p++)v=v+l[i][p]*r[p][k];l[i][k]=(a[i][k]-v)/r[k][k]; printf("l[%d][%d]=%f\n",i,k,l[i][k]);} for(i=1;i<n+1;i++){ v=0;for(k=1;k<i;k++)v=v+l[i][k]*Y[k];Y[i]=b[i]-v;}for(i=n;i>0;i--){ v=0;for(k=i+1;k<n+1;k++)v=v+r[i][k]*X[k];X[i]=(Y[i]-v)/r[i][i];}printf("\n");for(i=1;i<n+1;i++)printf("x[%d]=%f\n",i,X[i]);printf("\n");return ;}。

数值分析程序代码

数值积分h=0.01;x=0:h:1;y=4./(1+x.^2);format longt=length(x); %数组长度z1=sum(y(1:(t-1)))*h %矩形公式2z2=sum(y(2:t))*h %矩形公式2z3=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:nfor j=1:3X(j)=(b(j)-A(j,[1:j-1,j+1:3])*X0([1:j-1,j+1:3]))/A(j,j); endt=[t;X'];X0=X;enddisp('方程组的解')X1disp('迭代10步的解')Xdisp('方程组每次迭代的解')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:nfor j=1:3X(j)=(b(j)-A(j,[1:j-1,j+1:3])*X([1:j-1,j+1:3]))/A(j,j);endt=[t,X];enddisp('方程组的解')X1disp('迭代10步的解')Xdisp('方程组每次迭代的解')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）。

数值分析所有代码

实验一：拉格朗日插值多项式命名（源程序.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）=。

数值分析实验— MATLAB实现

数值分析实验——MATLAB实现姓名sumnat学号2013326600000班级13级应用数学2班指导老师2016年1月一、插值：拉格朗日插值 (1)1、代码： (1)2、示例： (1)二、函数逼近：最佳平方逼近 (2)1、代码： (2)2、示例： (2)三、数值积分：非反常积分的Romberg算法 (3)1、代码： (3)2、示例： (4)四、数值微分：5点法 (5)1、代码： (5)2、示例： (6)五、常微分方程：四阶龙格库塔及Adams加速法 (6)1、代码：四阶龙格库塔 (6)2、示例： (7)3、代码：Adams加速法 (7)4、示例： (8)六、方程求根：Aitken 迭代 (8)1、代码： (8)2、示例： (9)七、线性方程组直接法：三角分解 (9)1、代码： (9)2、示例： (10)八、线性方程组迭代法：Jacobi法及G-S法 (11)1、代码：Jacobi法 (11)2、示例： (12)3、代码：G-S法 (12)4、示例： (12)九、矩阵的特征值及特征向量：幂法 (13)1、代码： (13)2、示例： (13)一、插值：拉格朗日插值1、代码：function z=LGIP(x,y)%拉格朗日插值n=size(x);n=n(2);%计算点的个数syms a;u=0;%拉格朗日多项式f=1;%插值基函数for i=1:nfor j=1:nif j==if=f;elsef=f*(a-x(j))/(x(i)-x(j));endendu=u+y(i)*f;f=1;endz=expand(u);%展开2、示例：>> x=1:6;y1=x.^5+3*x.^2-6;y2=sin(x)+sqrt(x);>> f1=LGIP(x,y1)f1 =-6+3*a^2+a^5%可知多项式吻合得很好>> f2=vpa(LGIP(x,y2),3)f2 =.962e-1*a^4+1.38*a+.300*a^2+.504-.436*a^3-.616e-2*a^5二、函数逼近：最佳平方逼近1、代码：function z=BestF(u,a,b,n)%最佳平方逼近，用x^i逼近,n为逼近的次数n=n+1;syms xreal;old=findsym(u);u=subs(u,old,x); %将u中变量替换为xf=sym('');H=sym('');d=sym('');for i=1:n %生成函数系f(1,i)=x^(i-1);endfor i=1:n %生成内积Hfor j=1:nH(i,j)=int(f(1,i)*f(1,j),a,b);endendfor i=1:n %生成内积dd(i,1)=int(f(1,i)*u,a,b);enda=H\d;%解法方程Ha=dz=a'*f';2、示例：>> syms x>> f1=sqrt(x);>> f2=x^5+x^2;>> f3=exp(x);>> a=0 ;b=1;>> BestF(f1,a,b,5)ans =12/143+420/143*x-1120/143*x^2+2016/143*x^3-1800/143*x^4+56/13*x^5>> BestF(f2,a,b,5)ans =x^5+x^2>> BestF(f3,a,b,5)ans =-566826+208524*exp(1)+(16733010-6155730*exp(1))*x+(-115830120+42611520*exp(1))* x^2+(306348840-112699440*exp(1))*x^3+(-342469260+125987400*exp(1))*x^4+(136302012-5 0142708*exp(1))*x^5>> vpa(ans,3)ans =.1e4-.1e6*x-.1e7*x^3+.1e7*x^4三、数值积分：非反常积分的Romberg算法1、代码：function z=IntRom(f,a,b) %Romberg 算法e=1e-10;I{1}=linspace(a,b,2);%1等分I{2}=linspace(a,b,3);%2等分L=setdiff(I{2},I{1});%新得插值点h=b-a;T(1,1)=h/2*sum(subs(f,I{1}));T(2,1)=1/2*T(1,1)+h/2*sum(subs(f,L));T(2,2)=4/3*T(2,1)-1/3*T(1,1);k=2;while abs(T(k,k)-T(k-1,k-1))>e %精度要求k=k+1;I{k}=linspace(a,b,2^(k-1)+1);L=setdiff(I{k},I{k-1});%集合差运算，新得插值点h=h/2;T(k,1)=1/2*T(k-1,1)+h/2*sum(subs(f,L));%梯形for i=2:kT(k,i)=(4^(i-1)/(4^(i-1)-1))*T(k,i-1)-(1/(4^(i-1)-1))*T(k-1,i-1);%加速endEndz=T(k,k);2、示例：>> syms x>> f=x^4;>> a=-4;b=4;>> IntRom(f,a,b)T =1.0e+003 *2.04800000000000 0 0 01.02400000000000 0.68266666666667 0 00.57600000000000 0.42666666666667 0.40960000000000 00.45200000000000 0.41066666666667 0.40960000000000 0.40960000000000ans =4.096000000000000e+002>> vpa((int(f,a,b)-ans),3)ans =0.>> f=exp(x);>> a=0;b=1;>> IntRom(f,a,b)T =Columns 1 through 41.85914091422952 0 0 01.75393109246483 1.71886115187659 0 01.72722190455752 1.71831884192175 1.71828268792476 01.72051859216430 1.71828415469990 1.71828184221844 1.718281828794531.71884112857999 1.71828197405189 1.71828182867536 1.718281828460391.71842166031633 1.71828183756177 1.71828182846243 1.71828182845905Columns 5 through 60 00 00 00 01.71828182845908 01.71828182845905 1.71828182845905ans =1.71828182845905>> vpa((int(f,a,b)-ans),3)ans =0.四、数值微分：5点法1、代码：function z=VDiff(f,x0)%5点法求导数值e=1e-15;h=0.01;for i=0:4x(i+1)=x0+i*h;endy=subs(f,x);m(1)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));%5点导数公式h=h/2;for i=0:4x(i+1)=x0+i*h;endy=subs(f,x);m(2)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));h=h/2;for i=-0:4x(i+1)=x0+i*h;endy=subs(f,x);m(3)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));k=3;while abs(m(k)-m(k-1))<abs(m(k-1)-m(k-2)) & abs(m(k)-m(k-1))>e & (h/10)>0%控制收敛条件及精度要求及h非0h=h/2;k=k+1;for i=0:4x(i+1)=x0+i*h;endy=subs(f,x);m(k)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));ende=abs(m(k)-m(k-1));z=[m(k);e];2、示例：>> syms x>> f=exp(x);>> x0=2;>> VDiff(f,x0)ans =7.389056098949710.00000000002558五、常微分方程：四阶龙格库塔及Adams加速法1、代码：四阶龙格库塔function z=RGFour(f,y0,a,b)%4阶龙格库塔，f为函数句柄h=0.01;X=a:h:b;Y(1)=y0;n=size(X);n=n(2);for i=1:n-1K1=f([X(i) Y(i)]);K2=f([X(i)+h/2,Y(i)+h/2*K1]);K3=f([X(i)+h/2,Y(i)+h/2*K2]);K4=f([X(i)+h,Y(i)+h*K3]);Y(i+1)=Y(i)+h/6*(K1+2*K2+2*K3+K4);endz=Y;plot(X,Y);2、示例：>> f=@(x)sin(x(1));>> y0=0;a=0;b=2*pi;>> figure(1);>> RGFour(f,y0,a,b);3、代码：Adams加速法function z=myAdams(f,y0,a,b)h=0.01;p(4)=1;c(4)=1;X=a:h:b;Y(1)=y0;n=size(X);n=n(2);for i=1:3K1=f([X(i) Y(i)]);K2=f([X(i)+h/2,Y(i)+h/2*K1]);K3=f([X(i)+h/2,Y(i)+h/2*K2]);K4=f([X(i)+h,Y(i)+h*K3]);Y(i+1)=Y(i)+h/6*(K1+2*K2+2*K3+K4);endfor i=4:n-1p(i+1)=Y(i)+h/24*(55*f([X(i),Y(i)])-59*f([X(i-1),Y(i-1)])+37*f([X(i-2 ),Y(i-2)])-9*f([X(i-3),Y(i-3)]));m(i+1)=p(i+1)+251/720*(c(i)-p(i));m(i+1)=f([X(i+1),m(i+1)]);c(i+1)=Y(i)+h/24*(9*f([X(i+1),m(i+1)])+19*f([X(i),Y(i)])-5*f([X(i-1), Y(i-1)])+f([X(i-2),Y(i-2)]));Y(i+1)=c(i+1)-19/720*(c(i+1)-p(i+1));endz=Y;plot(X,Y);4、示例：>> f=@(x)exp(x(1));>> myAdams(f,0,0,2*pi);六、方程求根：Aitken 迭代1、代码：function z=myAitken(f,x0);%Aitken 迭代求方程的根e=1e-15;xu1=x0;xv1=subs(f,xu1);xv2=subs(f,xv1);if xv2-2*xv1+xu1==0%防止除0；xu2=xv2;elsexu2=xv2-(xv2-xv1)^2/(xv2-2*xv1+xu1);endwhile abs(xu2-xu1)>e%精度控制xu1=xu2;xv1=subs(f,xu1);xv2=subs(f,xv1);if xv2-2*xv1+xu1==0%防止除0；xu2=xv2;elsexu2=xv2-(xv2-xv1)^2/(xv2-2*xv1+xu1);%Aitken加速公式endendz=xu2;2、示例：>> syms x>> f=cos(x/2)+x;>> x0=3;>> myAitken(f,x0)ans =3.14159265358979>> f=x^2-2+x;>> x0=1;>> myAitken(f,x0)ans =1.41421356237309七、线性方程组直接法：三角分解1、代码：function z=myGuess(A,b);%线性方程组三角分解求根n=size(A);if n~=rank(A)z=['矩阵A线性相关，不符合要求'];return;endn=n(2);L=eye(n);for i=1:n-1for j=i+1:nL(j,i)=A(j,i)/A(i,i);A(j,:)=A(j,:)-L(j,i)*A(i,:);endendU=A;for i=1:n %解Ly=b，得ys=0;for j=1:i-1s=s+y(j)*L(i,j);endy(i)=(b(i)-s)/L(i,i);endfor i=n:-1:1 %解Ux=y，得xs=0;for j=i+1:ns=s+x(j)*U(i,j);endx(i)=(y(i)-s)/U(i,i);endLUz=x';2、示例：>> A=[1 2 3;2 1 5;11 17 21];>> b=[1 3 5]';>> myGuess(A,b)L =1.00000000000000 0 02.00000000000000 1.00000000000000 011.00000000000000 1.66666666666667 1.00000000000000U =1.000000000000002.000000000000003.000000000000000 -3.00000000000000 -1.000000000000000 0 -10.33333333333333ans =-0.06451612903226-0.580645161290320.74193548387097>> t=A\bt =-0.06451612903226-0.580645161290320.74193548387097八、线性方程组迭代法：Jacobi法及G-S法1、代码：Jacobi法function z=myJacobi(A,b)n=size(A);n=n(2);x{1}=zeros(n,1);%初始值e=1e-10;D=diag(diag(A));L=D-tril(A);U=D-triu(A);B=D\(L+U);f=D\b;Q=B'*B;[w,d]=eig(Q);p=max(abs(diag(d))) ; %谱半径if p>=1z=['迭代发散'];return;endx{2}=B*x{1}+f;k=2;while norm(x{k}-x{k-1})>ek=k+1;x{k}=B*x{k-1}+f;endz=x{k};2、示例：>> A=[8 -3 2;4 11 -1;6 3 12];>> b=[20 33 36]';>> myJacobi(A,b)ans =3.000000000013402.000000000012610.999999999992373、代码：G-S法function z=myGS(A,b)n=size(A);n=n(2);x{1}=zeros(n,1);e=1e-10;D=diag(diag(A));L=D-tril(A);U=D-triu(A);B=(D-L)\U;f=(D-L)\b;Q=B'*B;[w,d]=eig(Q);p=max(abs(diag(d))) ; %谱半径if p>=1z=['迭代发散'];return;endx{2}=B*x{1}+f;k=2;while norm(x{k}-x{k-1})>ek=k+1;x{k}=B*x{k-1}+f;endz=x{k};4、示例：>> A=[8 -3 2;4 11 -1;6 3 12];>> b=[20 33 36]';>> myGS(A,b)ans =3.000000000001351.999999999999160.99999999999954九、矩阵的特征值及特征向量：幂法1、代码：function z=myChar(A);%幂法求主特征值及对应的特征向量e=1e-10;n=size(A);n=n(2);v1=ones(n,1);u1=v1;v2=A*u1;a=min(v2);b=max(v2);if abs(a)>abs(b)c=a;elsec=b;endu2=v2/c;%规范化while norm(u2-u1)>eu1=u2;v2=A*u1;a=min(v2);b=max(v2);if abs(a)>abs(b)c=a;elsec=b;endu2=v2/c;%规范化endz{1}=c;z{2}=v2;2、示例：>> A=[8 -3 2;4 11 -1;6 3 12];>> u=myChar(A)u =[14.00000000046956] [3x1 double]>> u{1}ans =14.00000000046956 >> u{2}ans =4.20000000191478 0.93333333198946 14.00000000046956。

数值分析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)。

数值分析各种代码

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

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。