《常用数值分析算法.docx》由会员分享，可在线阅读，更多相关《常用数值分析算法.docx（47页珍藏版）》请在三一办公上搜索。
1、常用数值分析算法雅可比迭代求解方程组#include iostream.hvoid main int i,j,n;float *a,*x,e;void jacobi(float *a,float *x,float e,int n);cout<<请输入方程次数和要求精度:;cin>>n>>e; a=new float*n;x=new floatn;for(i=0;i<n;i+)ai=new floatn+2;cout<<请输入增广矩阵:<<endl;for(i=0;i<n;i+)for(j=0;j<n+1;j+)cin>>aij;cou
2、t<<请输入迭代初值:<<endl;for(i=0;i<n;i+)cin>>ain+1;jacobi(a,x,e,n); cout<<方程的解为：<<endl;for(i=0;i<n;i+)cout<<X<<i+1<<=<<xi<<endl; void jacobi(float *a,float *x,float e,int n)int i,j,flag;float k,t;for(i=0;i<n-1;i+) /选主元素for(flag=i,j=i+1;j<n;j+)if(aji>aflagi)
3、flag=j;if(flag!=i)for(j=0;j<n+1;j+)t=aflagj;aflagj=aij;aij=t;for(i=0;i<n;i+) /建立迭代矩阵for(k=aii,j=0;j<n+1;j+)aij/=k;while(flag) /迭代过程flag=0;for(i=0;i<n;i+)for(xi=ain,j=0;j<n;j+)if(j!=i)xi-=aij*ajn+1;if(xi-ain+1>e)flag=1;for(i=0;i<n;i+)ain+1=xi; 追赶法求解三对角方程组#include iostream.hvoid main int
4、 i,j,n;float *a;void catchwith(float *a,int n);cout<<请输入方程次数:;cin>>n; a=new float*n;for(i=0;i<n;i+)ai=new float5;cout<<请输入增广矩阵(只含主系数):<<endl;for(i=0;i<n;i+)for(j=0;j<4;j+)cin>>aij;catchwith(a,n); cout<<方程的解为：<<endl;for(i=0;i<n;i+)cout<<X<<i+1<<=<<a
5、i4<<endl; void catchwith(float *a,int n) for(int i;i<n;i+) /追ai3-=ai-13*ai0/ai-11;ai1-=ai-12*ai0/ai-11;ai0=0;an-14=an-13/an-11;for(i=n-2;i>=0;i+) /赶ai4=(ai3-ai2*ai+14)/an-11; 三次样条插值#include iostream.hvoid main int i,n;float *a,x;float three(float *a,int n,float x);cout<<请输入点的个数:;cin>&g
6、tn; a=new float*n;for(i=0;i<n;i+)ai=new float2; cout<<请输入各个点的坐标:<<endl;for(i=0;i<n;i+)cin>>ai0>>ai1;cout<<输入要求的自变量值:;cin>>x;cout<<方程的解为：<<three(a,n,x)<<endl; float three(float *a,int n,float x)void catchwith(float *a,int n);int i;float *b,y=0; /b保存三对角阵b=new fl
7、oat*n-1;for(i=0;i<n-1;i+)bi=new float5;for(i=0;i<n-1;i+) /初始化三对角阵bi0=1-(ai0-ai-10)/(ai+10-ai-10);bi1=2;bi2=1-bi0;bi3=3*(bi0*(ai1-ai-11)/(ai0-ai-10)+bi3*(ai+11-ai1)/(ai+10-ai0);b00=2,b01=b02,b02=0,bn-22=0;b03=b03-(1-b01)*(a11-a01)/(a10-a00);bn-23=bn-23-(1-bn-20)*(an-11-an-21)/(an-10-an-20);catchw
8、ith(b,n-1); /求解三对角阵for(i=0;x>ai0;i+); /寻找x的位置i-;y+=(x-ai0)/(ai+10-ai0)-1)*(x-ai0)/(ai+10-ai0)-1)*(2*(x-ai0)/(ai+10-ai0)+1)*ai1;y+=(x-ai0)/(ai+10-ai0)*(x-ai0)/(ai+10-ai0)*(-2*(x-ai0)/(ai+10-ai0)+3)*ai+11;y+=(ai+10-ai0)*(x-ai0)/(ai+10-ai0)*(x-ai0)/(ai+10-ai0)-1)*(x-ai0)/(ai+10-ai0)-1)*bi4;y+=(ai+10-
9、ai0)*(x-ai0)/(ai+10-ai0)*(x-ai0)/(ai+10-ai0)*(x-ai0)/(ai+10-ai0)-1)*bi+14;return y;void catchwith(float *a,int n) /追赶法解三对角阵 for(int i;i<n;i+) /追ai3-=ai-13*ai0/ai-11;ai1-=ai-12*ai0/ai-11;ai0=0;an-14=an-13/an-11;for(i=n-2;i>=0;i+) /赶ai4=(ai3-ai2*ai+14)/an-11; 高斯_塞德尔迭代求解方程组#include iostream.hvoid m
10、ain int i,j,n;float *a,e;void gaussSeidel(float *a,float e,int n);cout<<请输入方程次数和要求精度:;cin>>n>>e; a=new float*n;for(i=0;i<n;i+)ai=new floatn+2;cout<<请输入增广矩阵:<<endl;for(i=0;i<n;i+)for(j=0;j<n+1;j+)cin>>aij;cout<<请输入迭代初值:<<endl;for(i=0;i<n;i+)cin>>ain+1;gauss
11、Seidel(a,e,n); cout<<方程的解为：<<endl;for(i=0;i<n;i+)cout<<X<<i+1<<=<<ain+1<<endl; void gaussSeidel(float *a,float e,int n)int i,j,flag;float k,t;for(i=0;i<n-1;i+) /选主元素for(flag=i,j=i+1;j<n;j+)if(aji>aflagi)flag=j;if(flag!=i)for(j=0;j<n+1;j+)t=aflagj;aflagj=aij;aij=
12、t;for(i=0;i<n;i+) /建立迭代矩阵for(k=aii,j=0;j<n+1;j+)aij/=k;while(flag) /迭代过程flag=0;for(i=0;i<n;i+)for(k=ain,j=0;j<n;j+)if(j!=i)k-=aij*ajn+1;if(k-ain+1>e)flag=1;ain+1=k; 龙贝格求积#include <iostream>using namespace std;void mainfloat a,b,e,result; float f(float x);float romberg(float (*f)(float)
13、,float a,float b,float e);cout<<输入端点：;cin>>a>>b;cout<<输入精度：;cin>>e;result=romberg(f,a,b,e);cout<<运算结果是：<<result; float romberg(float (*f)(float),float a,float b,float e)/适合于比较光滑的曲线求积float t1,t2,s1,s2,c1,c2,r1=0,r2=e;float s,h,x,k=1;h=b-a,t1=h*(f(a)+f(b)/2;while(r2-r1>=e
14、 | r1-r2>=e)for(s=0,x=a+h/2;x<b;x+=h)s+=f(x);t2=(t1+h*s)/2;s2=t2+(t2-t1)/3;if(k=1)k+,h/=2;t1=t2,s1=s2;continue;c2=s2+(s2-s1)/15;if(k=2)k+,h/=2;t1=t2,s1=s2,c1=c2;continue;r2=c2+(c2-c1)/63;r1=r2,c1=c2;t1=t2,s1=s2;k+,h/=2;return r2;float f(float x)return x*x;中点加速求微#include <iostream>using names
15、pace std;void mainfloat x,e,result; float f(float x);float romberg(float (*f)(float),float x,float e);cout<<输入自变量值：;cin>>x;cout<<输入精度：;cin>>e;result=romberg(f,x,e);cout<<运算结果是：<<result<<endl; float romberg(float (*f)(float),float x,float e)/适合于比较光滑的曲线求积float t1,t2,h=1.0;f
16、loat g3(float x,float h);for(t2=g3(x,h),t1=t2-e;t1-t2>=e | t2-t1>=e;h/=2)t1=t2,t2=g3(x,h);return t2;/*/ /加速公式float g3(float x,float h)float g2(float x,float h);return (64*g2(x,h/2)-g2(x,h)/63;float g2(float x,float h)float g1(float x,float h);return (16*g1(x,h/2)-g1(x,h)/15;float g1(float x,floa
17、t h)float g(float x,float h);return (4*g(x,h/2)-g(x,h)/3;float g(float x,float h)float f(float x);return (f(x+h)-f(x-h)/2/h;/*/float f(float x)return x*x;变步长复化一点高斯求积#include <iostream>using namespace std;void mainfloat a,b,e,result; float f(float x);float guass(float (*f)(float),float a,float b,
18、float e);cout<<输入端点：;cin>>a>>b;cout<<输入精度：;cin>>e;result=guass(f,a,b,e);cout<<运算结果是：<<result; float guass(float (*f)(float),float a,float b,float e) /适合于比较光滑的曲线求积float t1,t2,s,h,x;h=(b-a)*2,t2=h*(f(a)+f(b)/2;for(t1=t2-2*e;t2-t1>=e | t1-t2>=e;t2=h*s) for(h/=2,t1=t2,s=0
19、,x=a+h/2;x<b;x+=h)s+=f(x); return t2;float f(float x)return x*x;复化辛甫生求积#include <iostream>using namespace std;void mainfloat a,b,result; int n;float f(float x);float simpson(float (*f)(float x),float a,float b,int n);cout<<输入端点：;cin>>a>>b;cout<<输入等分数：;cin>>n;result=simpson(
20、f,a,b,n);cout<<运算结果是：<<result; float simpson(float (*f)(float x),float a,float b,int n)float s,x,h;for(s=f(a)-f(b),x=a,h=(b-a)/n;x<b;x+=h)s+=4*f(x+h/2)+2*f(x+h);return h*s/6;float f(float x)return 17*x*x*x+79*x*x+23*x+107;快速弦截法求解方程#include <iostream>using namespace std;void mainfloat x
21、0,result,e;float f(float x);float f1(float x);float quick(float (*f)(float),float x0,float e);cout<<输入初始值：;cin>>x0;cout<<输入精度：;cin>>e;result=quick(f,x0,e);cout<<运算结果是：<<result<<endl;float quick(float (*f)(float),float x0,float e)float x1=x0+f(x0),x2;x2=x1-f(x1)*(x1-x0)/(f
22、(x1)-f(x0);while(x1-x2>=e | x2-x1>=e)x0=x1;x1=x2;x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0);return x2;float f(float x) /要解的方程为f(x)=0return x*x-x-12;亚当姆斯预报_校正系统解一阶微分的初值问题#include <iostream>using namespace std;void mainfloat x0,y0,h,*result;int N;float f(float x,float y);void adams(float (*f)(float,floa
23、t),float,float,int,float,float *result);cout<<输入初始点：;cin>>x0>>y0;cout<<输入解的个数和步长：;cin>>N>>h;result=new float *N;for(int i=0;i<N;i+)resulti=new float2;adams(f,x0,y0,N,h,result);cout<<运算结果是：<<endl; for(i=0;i<N;i+)cout<<resulti0<< <<resulti1<<endl;voi
24、d adams(float (*f)(float,float),float x0,float y0,int N,float h,float *result)float k1,k2,k3,k4;float yp,yp1;result00=x0;result01=y0;for(int n=1;n<N & n<4;x0+=h,y0=resultn+1) /经典格式求前三点resultn0=x0+h;k1=f(x0,y0);k2=f(x0+h/2,y0+h*k1/2);k3=f(x0+h/2,y0+h*k2/2);k4=f(x0+h,y0+h*k3);resultn1=y0+h*(k1+2*k
25、2+2*k3+k4)/6;for(;n<N;x0+=h,y0=resultn+1)resultn0=x0+h;yp=resultn-11+h*(55*f(resultn-10,resultn-11)-59*f(resultn-20,resultn-21)+37*f(resultn-30,resultn-31)-9*f(resultn-40,resultn-41)/24;yp1=f(resultn0,yp);resultn1=resultn-11+h*(9*yp1+19*f(resultn-10,resultn-11)-5*f(resultn-20,resultn-21)+f(resultn
26、-30,resultn-31)/24;float f(float x,float y) /当初始值满足y=x*x时，微分方程解为y=x*x;return 2*x+y-x*x;变步长复化梯形法求积#include <iostream>using namespace std;void mainfloat a,b,e,result; float f(float x);float exfloor(float (*f)(float),float a,float b,float e);cout<<输入端点：;cin>>a>>b;cout<<输入精度：;cin>&
27、gte;result=exfloor(f,a,b,e);cout<<运算结果是：<<result; float exfloor(float (*f)(float),float a,float b,float e) /对曲线光滑度要求不高float t1,t2,s,h,x;h=(b-a)*2,t2=h*(f(a)+f(b)/2;for(t1=t2-2*e;t2-t1>=e | t1-t2>=e;t2=(t1+h*s)/2) for(h/=2,t1=t2,s=0,x=a+h/2;x<b;x+=h)s+=f(x); return t2;float f(float x)re
28、turn x*x;改进的欧拉格式求解一阶微分的初值问题#include <iostream>using namespace std;void mainfloat x0,y0,h,*result;int N;float f(float x,float y);void euler(float (*f)(float,float),float,float,int,float,float *result);cout<<输入初始点：;cin>>x0>>y0;cout<<输入解的个数和步长：;cin>>N>>h;result=new float *N;
29、for(int i=0;i<N;i+)resulti=new float2;euler(f,x0,y0,N,h,result);cout<<运算结果是：<<endl; for(i=0;i<N;i+)cout<<resulti0<< <<resulti1<<endl;void euler(float (*f)(float,float),float x0,float y0,int N,float h,float *result)result00=x0;result01=y0;for(int n=1;n<N;x0+=h,y0=resultn
30、+1)resultn0=x0+h;resultn1=y0+h*(f(x0,y0)+f(x0+h,y0+h*f(x0,y0)/2;float f(float x,float y) /当初始值满足y=x*x时，微分方程解为y=x*x;return 2*x+y-x*x;经典格式式求解一阶微分的初值问题#include <iostream>using namespace std;void mainfloat x0,y0,h,*result;int N;float f(float x,float y);void runge_kutta(float (*f)(float,float),float,
31、float,int,float,float *result);cout<<输入初始点：;cin>>x0>>y0;cout<<输入解的个数和步长：;cin>>N>>h;result=new float *N;for(int i=0;i<N;i+)resulti=new float2;runge_kutta(f,x0,y0,N,h,result);cout<<运算结果是：<<endl; for(i=0;i<N;i+)cout<<resulti0<< <<resulti1<<endl;void r
32、unge_kutta(float (*f)(float,float),float x0,float y0,int N,float h,float *result)float k1,k2,k3,k4;result00=x0;result01=y0;for(int n=1;n<N;x0+=h,y0=resultn+1)resultn0=x0+h;k1=f(x0,y0);k2=f(x0+h/2,y0+h*k1/2);k3=f(x0+h/2,y0+h*k2/2);k4=f(x0+h,y0+h*k3);resultn1=y0+h*(k1+2*k2+2*k3+k4)/6;float f(float x
33、,float y) /当初始值满足y=x*x时，微分方程解为y=x*x;return 2*x+y-x*x;埃特金加速迭代求解方程#include <iostream>using namespace std;void mainfloat x0,result,e;float u(float x);float aitken(float (*f)(float),float x0,float e);cout<<输入初始值：;cin>>x0;cout<<输入精度：;cin>>e;result=aitken(u,x0,e);cout<<运算结果是：<<
34、result<<endl;float aitken(float (*u)(float),float x0,float e)float x1,x2;x1=u(x0),x2=u(x1);x2=x2-(x2-x1)*(x2-x1)/(x2-2*x1+x0);while(x2-x0>=e | x0-x2>=e)x0=x2;x1=u(x0);x2=u(x1);x2=x2-(x2-x1)*(x2-x1)/(x2-2*x1+x0);return x2;float u(float x) /要解的方程为x=u(x)return x*x-12;埃特金插值#include <iostream&g
35、tusing namespace std;void mainfloat *a,x,result,e; /保存精度int n,i;float aitken(float *a,int n,float x,float *e);cout<<输入点数：;cin>>n;a=new float *n;for(i=0;i<n;i+)ai=new float2;cout<<输入各个点：<<endl;for(i=0;i<n;i+)for(int j=0;j<2;j+)cin>>aij;cout<<输入自变量的值：;cin>>x;result=ai
36、tken(a,n,x,&e);cout<<运算结果是：<<result;if(e>0)cout<<+<<e<<endl;elsecout<<e<<endl;float aitken(float *a,int n,float x,float *e)for(int i=0;i<n-1;i+)for(int j=i+1;j<n;j+)aj1=ai1+(x-ai0)*(aj1-ai1)/(aj0-ai0);*e=an-11-an-21;return an-11;行列式求值#include <iostream>using nam
37、espace std;void mainfloat *a,result;int n,i;float hls(float *a,int n);cout<<输入行列式的阶数：;cin>>n;a=new float *n;for(i=0;i<n;i+)ai=new floatn;cout<<输入初始值：<<endl;for(i=0;i<n;i+)for(int j=0;j<n;j+)cin>>aij;result=hls(a,n);cout<<运算结果是：<<result<<endl;float hls(float *a,i
38、nt n)int i,j,k;float result=1;for(i=0;i<n-1;i+)for(j=i+1;j<n;j+)for(k=n;k>=0;k-)ajk=ajk-aik*aji/aii; for(i=0;i<n;i+)result=result*aii;return result;约当消去法求解方程组#include iostream.hvoid main int i,j,n;float *a;void jordan(float *a,int n);cout<<请输入方程次数:;cin>>n; a=new float*n;for(i=0;i<n
39、;i+)ai=new floatn+1;cout<<请输入增广矩阵:<<endl;for(i=0;i<n;i+)for(j=0;j<n+1;j+)cin>>aij;jordan(a,n); cout<<方程的解为：<<endl;for(i=0;i<n;i+)cout<<X<<i+1<<=<<ain<<endl; void jordan(float *a,int n)int i,j,k;for(i=0;i<n;i+)for(j=n;j>=i;j-)aij=aij/aii;for(k=0;k<n
40、;k+)if(k!=i)for(j=n;j>=i;j-)akj=akj-aij*aki; 最小二乘法#include <iostream>using namespace std;void mainfloat *a,x,result;int n;float k,b;void smallest(float *a,int n,float *k,float *b);cout<<输入点数：;cin>>n;a=new float *n;for(int i=0;i<n;i+)ai=new float2;cout<<输入各个点：<<endl;for(i=0;i<n;i+)for(in
