数值分析计算方法程序汇总.docx
- 文档编号:7297196
- 上传时间:2023-01-22
- 格式:DOCX
- 页数:21
- 大小:89.45KB
数值分析计算方法程序汇总.docx
《数值分析计算方法程序汇总.docx》由会员分享,可在线阅读,更多相关《数值分析计算方法程序汇总.docx(21页珍藏版)》请在冰豆网上搜索。
数值分析计算方法程序汇总
(一)秦九韶法
例:
已知p=2x2-3x-4,求x=1时,p=?
程序:
#include"stdio.h"
voidmain()
{
floata[50],b,x;
intn,i,k;
scanf("%d%f",&n,&x);
for(i=0;i<=n;i++)
scanf("%f",&a[i]);
b=a[0];
k=1;
while(k<=n)
{
b=x*b+a[k];
k=k+1;
}
printf("%f",b);
}
结果:
-5.000000
(二)复化中矩形公式
例:
求
=[g(y1)+g(y2)+···+g(ym)]h的值(已知g(y)=yn/2-1e-y/2,n=11,h=0.1,ym=20)。
程序:
#include"stdio.h"
#include
main()
{
doubley,h,s;
intn,i;
n=11;h=0.1;
s=0;y=0;
for(i=0;i<=200;i++)
{
y=y+h;
s=s+pow(y,n/2.0-1)*exp(-y/2.0);
}
s=s*h;
printf("%f\n",s);
}
结果:
2266.139467
(三)复化梯形公式
例:
求
程序:
#include"stdio.h"
voidmain()
{
doublea,b,s,h,x;
inti,n;
a=-1.0;b=1.0;n=10;
h=(b-a)/n;x=a;s=x*x/2;
for(i=1;i { x=x+h; s=s+x*x; } s=s+b*b/2; s=s*h; printf("s=%f\n",s); } 结果: s=0.680000 (四)复化辛普森公式 例: 求 程序: #include"stdio.h" voidmain() { doublea,b,c,s,h,x,y; inti,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 (五)复化高斯公式 例: 求 程序: #include #include main() {doublea,b,h,s,x1,x2; inti,n; a=0;b=2;n=20;s=0; h=(b-a)/n; for(i=0;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 (六)二维中矩形公式 例: 求 程序: #include #include main() {doublea,b,c,d,x[50],y[50],hx,hy,s; inti,j,n,k; a=0.0;b=1.0;c=0.0;d=1.0; n=k=10; hx=(b-a)/n;hy=(d-c)/k; x[0]=a+hx/2;y[0]=c+hy/2;s=0.0; for(j=0;j { for(i=0;i { s=s+x[i]*x[i]+y[j]*y[j]; x[i+1]=x[i]+hx; } y[j+1]=y[j]+hy; } s=hx*hy*s; printf("s=%f\n",s); } 结果: s=0.665000 (七)迭代法 例: 求x=x2的解。 程序: #include"stdio.h" #include main() { doublex,xl,y,yl; inti,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; elsexl=x; } for(j=0;;j++) {y=sqrt(y); if(fabs(yl-y)<0.0001) break; elseyl=y; } printf("x=%f,y=%f\n",x,y); } 结果: x=0.000000,y=0.999915 (八)牛顿迭代法 y=f(x),求f(x*)=0。 y-y0=f'(x0)(x-x0) -y0=f'(x0)(x1-x0) x1=x0-f(x0)/f'(x0) ....... xn+1=xn-f(xn)/f'(xn) 例: 求f(x)=x-x2的零点(已知x0=0.1)。 程序: #include"stdio.h" #include main() { doublex,xl; inti; x=0.1;xl=x; for(i=0;;i++) {x=x-(x-x*x)/(1-2*x); if(fabs(xl-x)<0.0001) break; elsexl=x; } printf("x=%f",x); } 结果: x=-0.000000 (九)二分法 已知f(a)f(b)<0,则在[a,b]中至少有一点f(x*)=0。 令x1=(a+b)/2 如f(a)f(x1)<0,取a1=a,b1=x1;否则,取a1=x1,b1=b。 令x2=(a1+b1)/2 如f(a1)f(x2)<0,取a2=a1,b2=x2;否则,取a2=x2,b2=b1。 ...... 令xn+1=(an+bn)/2 如f(an)f(xn+1)<0,取an+1=an,bn+1=xn+1;否则,取an+1=xn+1,bn+1=bn。 例: 求f(x)=x-x2在[0.1,2.45]内的零点。 程序: #include"stdio.h" #include main() { doublea,al,b,bl,x,xl; inti; a=0.1;b=2.45; for(i=0;;i++) {x=(a+b)/2; if((a-a*a)*(x-x*x)<0) al=a,b=x; elsea=x,bl=b; if(fabs(xl-x)<0.0001) break; elsexl=x; } printf("x=%f",xl); } 结果: x=1.000040 (十)一元回归 已知(x1,y1),(x2,y2)···(xN,yN), y=a+bxyi=a+bxi+EiQ(a,b)= Ei=yi-a-bxi 总误差Q= 令 得 aN+bX-Y=0,aX+bX2-XY=0 (其中X= Y= X2= XY= )。 可得: 例: 已知(0,0),(1,2),(2,5);求a,b。 程序: #include"stdio.h" main() { floata,b,E,X,Y,Xl,Xy,Q,x[3]={0.0,1.0,2.0},y[3]={0.0,2.0,5.0}; inti,j,n=3; X=0;Y=0;Xl=0;Xy=0; for(i=0;i {X=X+x[i]; Y=Y+y[i]; Xl=Xl+x[i]*x[i]; Xy=Xy+x[i]*y[i]; a=(Xl*Y-X*Xy)/(n*Xl-X*X); b=(n*Xy-X*Y)/(n*Xl-X*X); } printf("a=%f,b=%f\n",a,b); Q=0;E=0; for(j=0;j<=3;j++) {Q=Q+E*E; E=y[j]-a-b*x[j]; } printf("Q(总误差)=%f\n",Q); } 结果: a=-0.166667,b=2.500000 Q(总误差)=0.166667 (十一)泰勒插值 例: 计算当x=2.104时,ex=pn(x)=? 程序: #include"stdio.h" #include voidmain() {doublex,l,p; inti,n; x=2.104;n=100; p=1;l=1; for(i=0;i<=n;i++) { l=l*x/(i+1); p=p+l; } printf("p=%f\n",p); } 结果: p=8.198900 (十二)拉格朗日插值 例: 已知(0,1),(3,2),(8,3),计算x=5时,求p1(x)=? 程序: #include"stdio.h" #include voidmain() {floatx[3]={0.0,3.0,8.0},y[3]={1.0,2.0,3.0},xl,p,l; intn,k,j; n=2; xl=5.0;p=0.0; for(k=0;k<=n;k++) {l=1; for(j=0;j<=n;j++) {if(j==k) continue; elsel=l*(xl-x[j])/(x[k]-x[j]); } p=p+l*y[k]; } printf("p=%f\n",p); } 结果: p=2.500000 (十三)牛顿插值 已知(x0,y0),(x1,y1),···(xn,yn);求pn=a0+a1x+a2x2+···anxn ··························· (1) ······ 将上式依次代入 (1)中 例: 已知(0,0),(1,1),(2,4);求p2(x)=? 程序: #include #include main() {doublex[3]={0,1,2},y[3]={0,1,4},g[3],gl,p,pl,X; inti,n=2;X=4; g[0]=y[0]; g[1]=(y[1]-y[0])/(x[1]-x[0]); gl=(y[2]-y[0])/(x[2]-x[0]); g[2]=(gl-g[1])/(x[2]-x[1]); pl=g[n]; for(i=n;i>=1;i--) {p=pl*(X-x[i-1])+g[i-1]; pl=p; } printf("p=%f\n",p); } 结果: p=16.000000 (十四)分段插值 ······ 例: 程序: #include"stdio.h" #include voidmain() {doublex[10],y[10],p[10],h,xl; inti,n; n=10;x[0]=0.0;x[9]=2.0;xl=2.0; h=(x[9]-x[0])/n; scanf("%d",&i); x[i]=h*i; y[i]=exp(x[i]); x[i+1]=h*(i+1); y[i+1]=exp(x[i+1]); p[i+1]=y[i]+(y[i+1]-y[i])/(x[i+1]-x[i])*(xl-x[i]); printf("p[i+1]=%f\n",p[i+1]); } 结果: 输入i=0时,p[i+1]=3.214028 (十五)解线性方程组 例: 解方程组 程序: #include"stdio.h" voidmain() {doublea[3][3]={{5.0,1.0,1.0},{1.0,4.0,-1.0},{1.0,-1.0,6.0}},b[3]={1.0,2.0,3.0},x[3]; inti,j,k,n; n=3; for(k=0;k {for(j=k+1;j {for(i=k+1;i {a[i][j]=a[i][j]-a[k][j]*a[i][k]/a[k][k]; b[i]=b[i]-b[k]*a[i][k]/a[k][k]; } } } x[n-1]=b[n-1]/a[n-1][n-1]; for(k=n-2;k>=0;k--) {x[k]=(b[k]-a[k][k]*x[k+1]-a[k][n]*x[n])/a[k][k]; printf("%f\n",x[k]); } printf("%f\n",x[n-1]); } 结果: -0.414921 0.414921 0.572816 (十六)列主元消去法 例: 解方程组 程序: #include #include #defineN100 floata[N][N+1]; voidmain() {inti,j,k,n; floatt,s=0; printf("输入矩阵阶数: ");scanf("%d",&n); printf("\n");printf("输入增广矩阵: \n"); for(i=0;i for(j=0;j scanf("%f",&a[i][j]); for(k=0;k {for(i=k+1;i if(fabs(a[i][k])>fabs(a[k][k])) for(j=k;j {t=a[k][j];a[k][j]=a[i][j];a[i][j]=t;} for(i=k+1;i {a[i][k]=a[i][k]/a[k][k]; for(j=k+1;j } } a[n-1][n]=a[n-1][n]/a[n-1][n-1]; for(k=n-2;k>=0;k--) {s=0; for(j=k+1;j } for(i=0;i printf("x[%d]=%f\n",i+1,a[i][n]); } 结果: x[1]=2.000000 x[2]=1.000000 x[3]=-0.000000 (十七)追赶法 例: 程序: #include"stdio.h" voidmain() {doublex[3],a[4]={0,1,1,1},b[4]={4,4,4,4},c[4]={1,1,1,0},d[4]={1,2,1,0},e[4],f[4]; intk,n=3; e[0]=d[0]/b[0];f[0]=c[0]/b[0]; for(k=1;k<=n-1;k++) {e[k]=(d[k]-a[k]*e[k-1])/(a[k]*f[k-1]+b[k]); f[k]=c[k]/(a[k]*f[k-1]+b[k]); } x[n]=(d[n]-a[n]*e[n-1])/(b[n]-a[n]*f[n-1]); printf("%f\n",x[n]); for(k=n-1;k>=0;k--) {x[k]=e[k]-f[k]*x[k+1]; printf("%f\n",x[k]); } } 结果: -0.036900 0.147601 0.377035 0.155741 (十八)高斯-赛德尔迭代法 例: 程序: #include"stdio.h" #include #defineN100 voidmain() {doublea[N][N]={{8,-3,2},{4,11,-1},{6,3,12}},b[N]={20,33,36},x[N]={0,0,0},xl,s; inti,j,k,n=3; for(k=1;k<=10;k++) {for(i=0;i {s=0.0; xl=x[i]; for(j=0;j {if(j==i) continue; elses=s+a[i][j]*x[j]; } x[i]=(b[i]-s)/a[i][i]; } } for(i=0;i printf("x[%d]=%f\n",i+1,x[i]); } 结果: x[1]=3.000000 x[2]=2.000000 x[3]=1.000000 (十九)一步欧拉公式 例: 已知 求x=1时的值。 程序: #include"stdio.h" voidmain() { doublex,y,h; inti,n=10; h=0.1;x=0;y=1; for(i=1;i<=n;i++) { y=y+h*(x+y); x=x+h; } printf("%f\n",y); } 结果: 3.187485 (二十)隐式欧拉公式: 例: 已知 求x=1时的值。 程序: #include"stdio.h" #include voidmain() { doublex[10],y[10],yl,h,e; inti,j,n=10; h=0.1;e=0.01; x[0]=0.0;y[0]=1.0; for(i=0;i { y[i+1]=y[i]+h*(x[i]+y[i]); x[i+1]=x[i]+h; for(j=0;;j++) { y[j+1]=y[j]+h*(x[j+1]+y[j+1]); yl=y[j+1]; if(fabs(yl-y[i+1]) break; elsey[j+1]=yl; } } printf("%f\n",yl); } 结果: 3.728783 (二十一)四阶龙格-库塔公式 例: 已知y'=x+y,y(0)=1;求y(0.2)=? 程序: #include"stdio.h" voidmain() {doublex,y,h,k1,k2,k3,k4; inti,n; h=0.1;n=2;x=0.0;y=1.0; for(i=0;i { k1=x+y; k2=x+y+h/2+h*k1/2; k3=x+y+h/2+h*k2/2; k4=x+y+h+h*k3; y=y+h*(k1+2*k2+2*k3+k4)/6; x=x+h; } printf("y(0.2)=%f\n",y); } 结果: y(0.2)=1.242805 补充例题: (二十二) 例: 已知y''=x+y+y',用一步欧拉公式求n=10时,y=? ,y'=? 程序: #include"stdio.h" voidmain() { doublex,y,z,h; inti,n; h=0.1;n=10; x=0;y=0;z=1; for(i=1;i<=n;i++) { y=y+h*z; z=z+h*(x+y+z); x=x+h; } printf("y=%f,y'=%f\n",y,z); } 结果: y=2.001568,y'=4.314725 (二十三) 例: 已知y''=x+y,y(0)=1,y (1)=0,求n=3时,y (1)=? y (2)=? 程序: #include"stdio.h" voidmain() { doublex[3],y[3],d[3],e[3],f[3],b,h; inti,k,n=3; h=1.0/n;b=-2.0-h*h; x[1]=h;d[1]=x[1]*h*h-1.0; for(i=2;i { x[i]=x[i-1]+h; d[i]=h*h*x[i]; } e[1]=d[1]/b;f[1]=-1/b; for(k=2;k<=n-1;k++) { e[k]=(d[k]-e[k-1])/(b+f[k-1]); f[k]=-1.0/(b+f[k-1]); } y[n-1]=e[n-1]; printf("y[%d]=%f",n-1,y[n-1]); for(k=n-2;k>=1;k--) { y[k]=e[k]+f[k]*y[k+1]; printf("y[%d]=%f",k,y[k]); } } 结果: y[2]=0.233333y[1]=0.566667
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 计算方法 程序 汇总