C语言解题经典程序(2)

totoerer · 发表于 2004-5-15 07:11:13

<

center" align=center>/*列主元消去法解线性方程组 GAUSS.C */<

22pt; mso-line-height-rule: exactly">#include<math.h><

22pt; mso-line-height-rule: exactly">#include<stdio.h> int gauss(double **a,int n) { int i,j,k,l; double c,t,eps=1.0e-15; for(k=0;k<n;k++) { c=0.0; for(i=k;i<n;i++) { t=fabs(a[k]); if(t>c) { c=t;l=i; } } if(c<=eps)return(0); if(l!=k) for(j=k;j<=n;j++) { t=a[k][j]; a[k][j]=a[l][j]; a[l][j]=t; } c=1.0/a[k][k]; for(j=k+1;j<=n;j++)a[k][j]*=c; for(i=k+1;i<n;i++) for(j=k+1;j<=n;j++) a[j]-=(a[k]*a[k][j]); } for(i=n-1;i>=0;i--) for(j=i+1;j<n;j++) a[n]-=(a[j]*a[j][n]); return(1);} main() { int i,j,n=3; double *pa[3]; double a[3][4]={12,-3,3,15,-18,3,-1,-15,1,1,1,6}; for(i=0;i<n;i++)pa=a; if(gauss(pa,n)) { for(i=0;i<n;i++) printf("\nx%d = %15.14lf",i,a[n]); } else printf("\nFiled !"); }

totoerer · 发表于 2004-5-15 07:13:14

<

center" align=center>龙贝格求积公式计算定积分 Lomberg1.C<

24pt; mso-line-height-rule: exactly">#include<math.h><

24pt; mso-line-height-rule: exactly">#include<stdio.h>double f(double x) { return(4.0/(1+x*x)); }double lomberg(double a,double b,double eps) { int k=1; double h1,h2,t1,t2,s1,s2,c1,c2,r1,r2,sam; double rs,rc,rr,x; rs=1.0/3.0;rc=1.0/15.0;rr=1.0/63.0; h1=b-a; h2=h1*0.5; t1=(f(a)+f(b))*h2; t2=t1*0.5+f(a+h2)*h2; s1=(t2*4.0-t1)*rs; h1=h2;h2*=0.5;t1=t2;x=a+h2; t2=t1*0.5+(f(x)+f(x+h1))*h2; s2=(t2*4.0-t1)*rs; c1=(s2*16.0-s1)*rc; h1=h2;h2*=0.5;t1=t2;s1=s2;x=a+h2; sam=0.0; for(x=a+h2;x<b;x+=h1)sam+=f(x); t2=t1*0.5+sam*h2; s2=(t2*4.0-t1)*rs; c2=(s2*16.0-s1)*rc; r1=(c2*64.0-c1)*rr; while(k<=20) { h1=h2;h2*=0.5;t1=t2;s1=s2;c1=c2;x=a+h2; sam=0.0; for(x=a+h2;x<b;x+=h1)sam+=f(x); t2=t1*0.5+sam*h2; s2=(t2*4.0-t1)*rs; c2=(s2*16.0-s1)*rc; r2=(c2*64.0-c1)*rr; if(fabs(r2-r1)<eps)break; r1=r2;k++; } printf("\nk=%d",k); return(r2); }main() { double eps=1.0e-12,answer; answer=lomberg(0.0,1.0,eps); printf("\nThe Answer is %15.14lf ",answer); }循环3次计算结果: 3.14159265358979

totoerer · 发表于 2004-5-15 07:13:40

<

center" align=center> 步长自适应龙格—库塔法解微分方程 RUNGKUTA.C<

22pt; mso-line-height-rule: exactly">#include<math.h><

22pt; mso-line-height-rule: exactly">#include<stdio.h> double f(double x,double y){ return(x*exp(-y));} double rungkuta(double x0,double y0,double h){double h12,h16,x,y,K1,K2,K3,K4;h12=0.5*h;h16=h/6.0; x=x0;y=y0;K1=f(x,y); x+=h12;y=y0+h12*K1;K2=f(x,y); y=y0+h12*K2;K3=f(x,y); x+=h12;y=y0+h*K3;K4=f(x,y); y=y0+(K1+2.0*K2+2.0*K3+K4)*h16;return(y);} double get_step(double x0,double y0,double *h,double eps){int k=0; double h1,y,y1,y2;h1=*h; y=rungkuta(x0,y0,h1);while(k<32) { h1*=0.5; y1=rungkuta(x0,y0,h1); y2=rungkuta(x0+h1,y1,h1); if(fabs(y2-y)<eps)break; y=y1;k++; }if(k==0)return(y2);*h=h1;return(y1);}double r_k(double x0,double y0,double xn,double eps){double h,x,y;h=(xn-x0)/256.0; x=x0;y=y0;for(;;) { y=get_step(x,y,&h,eps); x+=h; if(x>=xn)break; }return(y);} main(){double answer,eps=1.0e-15;printf("\nThe Answer is %16.15lf",log(13.5));answer=r_k(0.0,0.0,5.0,eps);printf("\nThe Answer is %16.15lf",answer);}计算结果：<v:shapetype> <v:stroke joinstyle="miter"></v:stroke><v:formulas><v:f eqn="if lineDrawn pixelLineWidth 0"></v:f><v:f eqn="sum @0 1 0"></v:f><v:f eqn="sum 0 0 @1"></v:f><v:f eqn="prod @2 1 2"></v:f><v:f eqn="prod @3 21600 pixelWidth"></v:f><v:f eqn="prod @3 21600 pixelHeight"></v:f><v:f eqn="sum @0 0 1"></v:f><v:f eqn="prod @6 1 2"></v:f><v:f eqn="prod @7 21600 pixelWidth"></v:f><v:f eqn="sum @8 21600 0"></v:f><v:f eqn="prod @7 21600 pixelHeight"></v:f><v:f eqn="sum @10 21600 0"></v:f></v:formulas><v:path gradientshapeok="t" connecttype="rect" extrusionok="f"></v:path><lock v:ext="edit" aspectratio="t"></lock></v:shapetype><v:shape><v:imagedata></v:imagedata></v:shape>精确结果为：<v:shape> <v:imagedata></v:imagedata></v:shape>

totoerer · 发表于 2004-5-15 07:14:18

<

center" align=center>幂法求矩阵的主特征值tezheng.c<

25pt; mso-line-height-rule: exactly">#include<stdio.h><

25pt; mso-line-height-rule: exactly">#include<math.h> double mifa(double *a[],double u0[],int n,double eps) { int k=1,i,j; double t0,t1,u1[20]; t0=u0[0]; for(i=1;i<n;i++) if(fabs(u0)>fabs(t0))t0=u0; while(k<=200) { for(i=0;i<n;i++) { u1=0.0; for(j=0;j<n;j++)u1+=a[j]*u0[j]; } t1=u1[0];for(i=1;i<n;i++)if(fabs(u1)>fabs(t1)) t1=u1; for(i=0;i<n;i++) u1/=t1; for(i=0;i<n;i++)u0=u1; if(fabs(t1-t0)<eps)break; k++;t0=t1; } return(t1); } main() { int i,n=3; double a[3][3]={6,-12,6,-21,-3,24,-12,-12,51}; double *pa[3],u0[3]={1,2,3}; double lbt,eps=1.0e-12; pa[0]=a[0];pa[1]=a[1];pa[2]=a[2]; lbt=mifa(pa,u0,n,eps); printf("\nlbt=%10.8lf",lbt); printf("\n"); for(i=0;i<n;i++)printf("X[%d]=%10.8lf ",i,u0); }

totoerer · 发表于 2004-5-15 07:15:07

<

center" align=center>牛顿插值法通用计算程序NUCA.C<

24pt; mso-line-height-rule: exactly">#include<stdio.h><

24pt; mso-line-height-rule: exactly">#include<conio.h>double nuca(double xi[],double yi[],int n,double x){int k,j; double y,xn; double yk[40];for(k=0;k<=n;k++)yk[k]=yi[k];y=yk[0]; xn=1.0;for(j=1;j<=n;j++) { for(k=0;k<=n-j;k++) yk[k]=(yk[k+1]-yk[k])/(xi[k+j]-xi[k]); xn=xn*(x-xi[j-1]); y+=yk[0]*xn; }return(y);}main(){double xi[]={64.0,81.0,100.0,121.0,144.0,169.0};double yi[]={8.0,9.0,10.0,11.0,12.0,13.0};double x＝115.0,y;y=nuca(xi,yi,5,x);printf("\nx=%g y=%10.8lf",x,y);}

totoerer · 发表于 2004-5-15 07:16:04

<

center" align=center>牛顿法解非线性方程组 newtonfa.c<

23pt; mso-line-height-rule: exactly">#include<math.h><

23pt; mso-line-height-rule: exactly">#include<stdio.h>#define N 3double a[N][N+1];void FDF(double x[]){double ex,sx,pi;pi=4.0*atan(1.0);a[0][3]=3.0*x[0]-cos(x[1]*x[2])-0.5;a[1][3]=x[0]*x[0]-81.0*(x[1]+0.1)*(x[1]+0.1)+sin(x[2])+1.06;ex=exp(-x[0]*x[1]);a[2][3]=ex+20.0*x[2]+pi*10.0/3.0-1.0;a[0][0]=3.0;sx=sin(x[1]*x[2]);a[0][1]=x[2]*sx;a[0][2]=x[1]*sx;a[1][0]=2.0*x[0];a[1][1]=-162.0*(x[1]+0.1);a[1][2]=cos(x[2]);a[2][0]=-x[1]*ex;a[2][1]=-x[0]*ex;a[2][2]=20.0;}int gauss(int n) { int i,j,k,i0; double c,t,eps=1.0e-15; for(k=0;k<n;k++) { c=0.0; for(i=k;i<n;i++) { t=fabs(a[k]); if(t>c) { c=t;i0=i; } } if(c<=eps)return(0); if(i0!=k) for(j=k;j<=n;j++) { t=a[k][j]; a[k][j]=a[i0][j]; a[i0][j]=t; } c=1.0/a[k][k]; for(j=k+1;j<=n;j++)a[k][j]*=c; for(i=k+1;i<n;i++) for(j=k+1;j<=n;j++) a[j]-=(a[k]*a[k][j]); } for(i=n-1;i>=0;i--) for(j=i+1;j<n;j++) a[n]-=(a[j]*a[j][n]); return(1);}int newtonfa(double x0[],double eps){int i,k=1; double dx,x1[N],y[N];for(;;) { FDF(x0); if(!gauss(N))return(0); for(i=0;i<N;i++)y=a[N]; dx=0.0; for(i=0;i<N;i++) if(fabs(y)>dx)dx=fabs(y); for(i=0;i<N;i++) x1=x0-y; for(i=0;i<N;i++) x0=x1; if(dx<eps)break; printf("\nk=%d",k); printf("\nx1=%12.10lf x2=%12.10lf x3=%12.10lf",x1[0],x1[1],x1[2]); k++;if(k>200)return(0); }return(1);}main(){int i;double x0[N]={3.0,2.0,1.0},eps=1.0e-12;if(newtonfa(x0,eps)) { for(i=0;i<N;i++) printf("\nx[%d] = %14.12lf ",i,x0); }else printf("\nFiled !");}

totoerer · 发表于 2004-5-15 07:17:08

<

center" align=center>牛顿法求方程的根 NTROOT.C<

22pt; mso-line-height-rule: exactly">#include<math.h><

22pt; mso-line-height-rule: exactly">#include<stdio.h>double f(double x){ return(x+log(x)-2.0); } double df(double x){ return(1.0+1.0/x); } double ca_root(double x0,double eps){int n=1; double x1,fx,dfx;while(n<800) { fx=f(x0); dfx=df(x0); x1=x0-fx/dfx; printf("\nn=%d x1=%14.12lf",n,x1); if(fabs(x1-x0)<eps)break; x0=x1;n++; }return(x1);}main(){double eps=1.0e-12,root;root=ca_root(1.0,eps);printf("\nThe root =%14.12lf",root);}经过5次循环，计算结果为：1.557145598998

totoerer · 发表于 2004-5-15 07:21:26

<

>/*穷举法求多元函数极值*/<

>#include<math.h><

>#include<float.h>#include<stdio.h>#include<stdlib.h>double PI,xishu;double ex=600.0,dx=190.0;double N01(double t){double x;x=(t-ex)/dx;return(exp(-(x*x*0.5)));}double EXN01(double t){double x;x=(t-ex)/dx;return(t*exp(-x*x*0.5));}double simpson(double a,double b,double eps,double (*f)(double)) { int m=1; double x,h,h1,sam1,sam2,s2,s1,d,rs; h=0.5*(b-a); rs=1.0/3.0; sam1=f(a)+f(b); sam2=f(a+h); s1=h*rs*(sam1+4.0*sam2); while(m<=19) { h1=h; h*=0.5; sam1+=(sam2*2.0); sam2=0.0; x=a+h; while(x<b) { sam2+=f(x); x+=h1; } s2=h*rs*(sam1+4.0*sam2); if(fabs(s2-s1)<eps)break; s1=s2; m++; } return(s2); }double F1(double x) { return(xishu*simpson(0.0,x,1.0e-12,EXN01)); }double FX(double x) { return(xishu*simpson(0.0,x,1.0e-12,N01)); }double au(double x) { double a,g1,g2; g1=F1(x); g2=FX(x); a=(g1+x-x*g2)/g2; return(a); } main() { double minl=1.0e+10,x,y,z,a,b=11400.0,c,f2; double minx,miny; PI=4.0*atan(1.0); xishu=1.0/(dx*sqrt(PI*2.0)); for(x=300.0;x<=400.0;x+=1.0) { f2=1000.0/x; a=au(x);c=(a*b)/(a+b); for(y=10.0;y<=25.0;y+=1.0) { z=12.0/y+f2+((y+1.0)*100.0+2500)/c; if(z<minl) { minl=z;minx=x,miny=y; printf("\nminl=%10.8lf x=%g y=%g",minl,minx,miny); } } } printf("\nminl=%10.8lf x=%g y=%g",minl,minx,miny); }

totoerer · 发表于 2004-5-15 07:22:16

<

36pt; TEXT-INDENT: -36pt; LINE-HEIGHT: 22pt; mso-line-height-rule: exactly; tab-stops: list 36.0pt; mso-list: l2 level1 lfo2">例3 当n=0,1,2,…,10时，输出3n<

22pt; mso-line-height-rule: exactly">#include<stdio.h><

22pt; mso-line-height-rule: exactly">main() { int n; unsigned answer=1;for(n=0;n<=10;n++) { printf(“\n3的%d次方等于%u”,n,answer); answer*=3;}} 例4 输入两个任意正整数，求它们的最大公约数。#include<stdio.h>main(){ unsigned n1, n2, t;printf(“\nInput two number\n”);scanf(“%u,%u”,&n1,&n2);if(n1>n2) {t=n1; n1=n2; n2=t;}for(;;){t=n2%n1;if(t==0)break;n2=n1; n1=t;}printf(“\n最大公约数是%u”,n1);}

totoerer · 发表于 2004-5-15 07:22:46

<

center" align=center>/* 列主元法求逆矩阵 */<

18pt">#include<math.h><

18pt">#include<stdio.h>#include<malloc.h>int inverta(double *a[],double *b[],int n,double eps) { int i,j,k,i0; double c,t,maxa; for(i=0;i<n;i++) for(j=0;j<n;j++) if(j==i)b[j]=1.0; else b[j]=0.0; for(k=0;k<n;k++){maxa=0.0; for(i=k;i<n;i++) { t=fabs(a[k]); if(t>maxa) { maxa=t;i0=i; } } if(maxa<=eps)return(0); if(i0!=k) { for(j=0;j<n;j++) { t=a[k][j];c=b[k][j]; a[k][j]=a[i0][j];b[k][j]=b[i0][j]; a[i0][j]=t;b[i0][j]=c; } } c=1.0/a[k][k]; for(j=k;j<n;j++)a[k][j]*=c; for(j=0;j<n;j++)b[k][j]*=c; for(i=k+1;i<n;i++) { for(j=k+1;j<n;j++) a[j]-=(a[k]*a[k][j]); for(j=0;j<n;j++) b[j]-=(a[k]*b[k][j]); } } for(k=n-2;k>=0;k--) for(i=k;i>=0;i--) for(j=0;j<n;j++)b[j]-=(b[k+1][j]*a[k+1]); return(1); } main() { int i,j,n=4; double eps=1.0e-9, *pa[4],*pb[4]; double a[4][4]={1,1,0,1,1,2,2,2,2,-2,2,1,3,-1,5,3},b[4][4]; for(i=0;i<n;i++) { pa=a; pb=b; } if(inverta(pa,pb,n,eps)) { for(i=0;i<n;i++) { printf("\n"); for(j=0;j<n;j++) printf("a-1[%d][%d]=%g ",i,j,b[j]); } } else printf("\nFiled !"); }

		自动登录	找回密码
密码			注-册-帐-号