本文精选了多种C语言实现的数值算法案例，包括显示方法、艾特肯加速收敛法、杜氏分解法等，覆盖二分法、插值、积分、微分方程求解等核心数值计算方法。
C数值算法
关键词： c程序 二分法 改进欧拉法
1、C程序
（1）I显示方法
#include <stdio.h>
#define HEIGHT 17
int main(void)
{
int i=0;
printf("/n/nIIIIIII/n");
while (i<HEIGHT)
{
printf(" III/n");
i++;
}
printf("IIIIIII/n/n/n");
return 0;
}
（2）艾特肯
#include<stdio.h>
#include<math.h>
float f(float a)
{float b;
b=pow(2,-a);
return b;
}
main()
{float x0,x1,e,y,z;
int k;
scanf("%f,%f",&x0,&e);
for(k=0;k<100;k++)
{y=f(x0);
z=f(y);
x1=x0-(y-x0)*(y-x0)/(z-2*y+x0);
if(fabs(x1-x0)<e)
{printf("%f",x1);
break;
}
else x0=x1;
}
}
（3）杜氏分解法
#include<stdio.h>
#include<math.h>
main()
{int n,i,j,k,t,sum;
float a[10][10],b[10],x[10],y[10];
printf("the top exp is ");
scanf("%d",&n);
for(i=0;i<n;i++)
for(j=0;j<n;j++)
scanf("%d",&a[i][j]);
for(i=0;i<n;i++)
scanf("%d",&b[i]);
for(i=1;i<n;i++)
a[i][1]=a[i][1]/a[1][1];
for(k=1;k<n;k++)
{for(j=k;j<n;j++)
{sum=0;
for(t=0;t<k;t++)
sum+=a[k][t]*a[t][j];
a[k][j]=a[k][j]-sum;
}
if(k!=n-1)
{for(i=k+1;i<n;i++)
{sum=0;
for(t=0;t<k;t++)
sum+=a[i][t]*a[t][k];
a[i][k]=(a[i][k]-sum)/a[k][k];
}
}
}
y[0]=b[0];
for(i=1;i<n;i++)
{sum=0;
for(k=0;k<i;k++)
sum+=a[i][k]*y[k];
y[i]=b[i]-sum;
}
x[n-1]=y[n-1]/a[n-1][n-1];
for(i=n-2;i>=0;i--)
{sum=0;
for(k=i+1;k<n;k++)
sum+=a[i][k]*x[k];
x[i]=(y[i]-sum)/a[i][i];
}
for(i=0;i<n;i++)
printf("%f/n",x[i]);
}
（4）二分法
#include<stdio.h>
#include<math.h>
float f(float x)
{float y;
y=x*x+2*x+1;
return y;
}
main()
{float a,b,c,x;
printf("%s","input three numbers:/n");
scanf("%f%f%f",&a,&b,&c);
while(fabs(b-a)>=c)
{x=(a+b)/2;
if(f(x)!=0)
{if(f(a)*f(x)<0)
b=x;
else a=x;}
else a=x,b=x;
}
printf("%f/n",(a+b)/2);
}
（5）分段线性插值
#include<stdio.h>
#include<math.h>
main()
{int e,h,i,j,k,n;
float x[10],y[10],a,sum,total,e1,e2;
printf("input n/n");
scanf("%d",&n);
printf("input X/n");
for(i=0;i<=n;i++)
scanf("%f",&x[i]);
printf("input Y/n");
for(i=0;i<=n;i++)
scanf("%f",&y[i]);
printf("input x/n");
scanf("%f",&a);
for(j=1;j<=n;j++)
if(x[j-1]<=a&&a<=x[j])
e=j;
if(e==1)
h=j;
else if(e==n)
h=j-1;
else {e1=fabs(a-x[j-2]);
e2=fabs(a-x[j+1]);
if(e1<e2)
h=j-1;
else h=j;}
sum=0;
for(k=h-1;k<=h+1;k++)
{total=1;
for(j=h-1;j<k;j++)
total*=((a-x[j])/(x[k]-x[j]));
for(j=k+1;j<=h+1;j++)
total*=((a-x[j])/(x[k]-x[j]));
sum+=y[k]*total;
}
printf("x=%f,L=%f",a,sum);
}
（6）复合梯形法
#include<stdio.h>
#include<math.h>
#define f(x) 4/(1+(x)*(x))
float g(int k)
{int i,total;
total=1;
for(i=1;i<=k;i++)
total*=2;
return total;
}
main()
{int i,k,k0,k1,tag;
float a,b,e,e1,t[50],sum;
scanf("%f,%f,%f,%d,%d",&a,&b,&e,&k0,&k1);
t[0]=((b-a)*(f(1)-f(0)))/2;
k=1;
tag=0;
do{sum=0;
for(i=1;i<=g(k-1);i++)
sum+=f(a+(2*i-1)*(b-a)/g(k));
t[k]=t[k-1]/2+(b-a)*sum/g(k);
if(k>=k1)
{e1=fabs(t[k]-t[k-1]);
if(e1<e)
{printf("k=%d,T=%f/n",k,t[k]);
tag=1;
}
}
k++;
}while(tag==0&&k<k0);
if(tag==0||k>=k0)
printf("error!/n");
}
（7）复合辛普森
#include<stdio.h>
#include<math.h>
#define f(x) (x)/(4+(x)*(x))
main()
{int k,n;
float a,b,h,s,x;
scanf("%f,%f,%d",&a,&b,&n);
h=(b-a)/(2*n);
x=a;
s=f(x)-f(b);
k=1;
do
{x=x+h;s=s+4*f(x);
x=x+h;s=s+2*f(x);
k++;
}while(k<=n);
s=(s*h)/3;
printf("s=%f/n",s);
}
（8）改进欧拉法
#include<stdio.h>
#include<math.h>
#define f(x,y) (x)-(y)+1
main()
{int k,n;
float a,b,h,y0,x,y,t[2];
scanf("%f,%f,%f,%d",&a,&b,&y0,&n);
h=(b-a)/n;
x=a;y=y0;
printf("%f,%f/n",x,y);
k=1;
do{t[0]=y+h*f(x,y);
x=x+h;
t[1]=y+h*f(x,t[0]);
y=(t[0]+t[1])/2;
printf("%f,%f/n",x,y);
k++;
}while(k<=n);
}
（9）高斯消去法
#include<stdio.h>
#include<math.h>
main()
{float a[10][10],b[10],m[10][10],x[10],sum;
int i,j,k,n;
printf("the top exp:");
scanf("%d",&n);
printf("/n");
for(i=0;i<n;i++)
for(j=0;j<n;j++)
scanf("%f",&a[i][j]);
for(i=0;i<n;i++)
scanf("%f",&b[i]);
for(k=0;k<n-1;k++)
{if(a[k][k]==0)
printf("error");
else for(i=k+1;i<n;i++)
{m[i][k]=a[i][k]/a[k][k];
a[i][k]=m[i][k];
b[i]=b[i]-m[i][k]*b[k];
for(j=k+1;j<n;j++)
a[i][j]=a[i][j]-m[i][k]*a[k][j];
}}
if(a[n-1][n-1]==0)
printf("error");
else x[n-1]=b[n-1]/a[n-1][n-1];
b[n-1]=x[n-1];
for(i=n-2;i>=0;i--)
{sum=0;
for(j=i+1;j<n;j++)
{sum+=a[i][j]*x[j];}
x[i]=(b[i]-sum)/a[i][i];
b[i]=x[i];
}
for(i=0;i<n;i++)
printf("%f/n",x[i]);
}
（10）简单迭代法
#include<stdio.h>
#include<math.h>
float f(float x)
{float y;
y=sqrt(10/(4+x));
return y;
}
main()
{float x0,x1,a;
int k,N;
k=0;
scanf("%f,%f,%d",&x0,&a,&N);
for(k=0;k<N;k++)
{x1=f(x0);
if(fabs(x1-x0)<a)
{printf("%d,%f",k,x1);
break;
}
else x0=x1;
}
if(k>=N)
{printf("error");
}
}
（11）列主元元素消元
#include<stdio.h>
#include<math.h>
main()
{float a[10][10],b[10],s,t,e,sum;
int i,j,k,n,m;
printf("The top exp is ");
scanf("%d",&n);
for(i=0;i<n;i++)
for(j=0;j<n;j++)
scanf("%f",&a[i][j]);
for(i=0;i<n;i++)
scanf("%f",&b[i]);
scanf("%f",&e);
k=0;
do{t=a[k][k];
for(i=k;i<n;i++)
{if(fabs(t)<fabs(a[i][k]))
{t=a[i][k];
m=i;
}
else m=k;
}
if(fabs(t)<e)
printf("det A = 0/n");
else {if(m!=k)
{for(j=0;j<n;j++)
{s=a[m][j];
a[m][j]=a[k][j];
a[k][j]=s;
}
s=b[m];
b[m]=b[k];
b[k]=s;
}
for(i=k+1;i<n;i++)
for(j=k+1;j<n;j++)
{a[i][k]=a[i][k]/a[k][k];
a[i][j]=a[i][j]-a[i][k]*a[k][j];
b[i]=b[i]-a[i][k]*b[k];
}
}
k++;
}while(k<n-2);
if(fabs(a[n-1][n-1])<e)
printf("det A = 0/n");
else {b[n-1]=b[n-1]/a[n-1][n-1];
for(i=n-2;i>=0;i--)
{sum=0;
for(k=i+1;k<n;k++)
{sum+=a[k][j]*b[j];}
b[i]=(b[i]-sum)/a[i][i];
}
}
for(i=0;i<n;i++)
printf("%f/n",b[i]);
}
（12）龙贝格算法
#include<stdio.h>
#include<math.h>
#define f(x) 4/(1+(x)*(x))
float g(int k)
{int i,total;
total=1;
for(i=1;i<=k;i++)
total*=2;
return total;
}
main()
{int k0,k1,k,tag,i;
float a,b,e,e1,t[50],s[50],c[50],r[50],sum;
scanf("%f,%f,%f,%d,%d",&a,&b,&e,&k0,&k1);
t[0]=(b-a)*(f(0)+f(1))/2;
k=1;
tag=0;
do{if(k==1)
{t[1]=t[0]/2+(b-a)*f(a+(b-a)/2)/2;
s[0]=(4*t[1]-t[0])/3;
}
if(k==2)
{t[2]=t[1]/2+(b-a)*(f(a+(b-a)/4)+f(a+3*(b-a)/4))/4;
s[1]=(4*t[2]-t[1])/3;
c[0]=(16*s[1]-s[0])/15;
}
if(k>=3)
{sum=0;
for(i=1;i<=g(k-1);i++)
sum+=f(a+(2*i-1)*(b-a)/g(k));
t[k]=t[k-1]/2+(b-a)*sum/g(k);
s[k-1]=(4*t[k]-t[k-1])/3;
c[k-2]=(16*s[k-1]-s[k-2])/15;
r[k-3]=(64*c[k-2]-c[k-3])/63;
}
if(k>=k1)
{e1=fabs(r[k-3]-r[k-4]);
if(e1<e)
{printf("k=%d,T=%f/n",k,r[k-3]);
tag=1;
}
}
k++;
}while(tag==0&&k<k0);
if(tag==0||k>=k0)
printf("error!/n");
}
（13）龙格库塔方法
#include<stdio.h>
#include<math.h>
#define f(x,y) (x)-(y)+1
main()
{int n,k;
float a,b,y0,h,x,y,t[4];
scanf("%f,%f,%f,%d",&a,&b,&y0,&n);
h=(b-a)/n;
x=a;y=y0;
printf("%f,%f/n",x,y);
k=1;
do{t[0]=f(x,y);
x=x+h/2;
t[1]=f(x,y+h*t[0]/2);
t[2]=f(x,y+h*t[1]/2);
x=x+h/2;
t[3]=f(x,y+h*t[2]);
y=y+h*(t[0]+2*t[1]+2*t[2]+t[3])/6;
printf("%f,%f/n",x,y);
k++;
}while(k<=n);
}
（14）牛顿插值多项式
#include<stdio.h>
#include<math.h>
main()
{int i,j,k,n;
float x[10],y[10],a,b,p;
printf("input n/n");
scanf("%d",&n);
printf("input X/n");
for(i=0;i<=n;i++)
scanf("%f",&x[i]);
printf("input Y/n");
for(i=0;i<=n;i++)
scanf("%f",&y[i]);
printf("input x/n");
scanf("%f",&a);
b=0;
k=0;
do{p=1;
j=0;
do
{if(j!=k)
p=p*(a-x[j])/(x[k]-x[j]);
j++;
}while(j<n);
b=b+p*y[k];
k++;
}while(k<n);
printf("x=%f,y=%f",a,b);
}
（15）牛顿迭代
#include<stdio.h>
#include<math.h>
#define N 100
float f(float x)
{float y;
y=x-cos(x);
return y;
}
main()
{float x0,x1,a;
int k;
scanf("%f,%f",&x0,&a);
for(k=0;k<N;k++)
{x1=x0-(x0-cos(x0))/(1+sin(x0));
if(fabs(x1-x0)<a)
{printf("%f,%d",x1,k);
break;
}
else x0=x1;
}
if(k>=N)
{printf("error");
}
}
（16）牛顿下山
#include<stdio.h>
#include<math.h>
float f(float x)
{float y;
y=pow(x,3)-x-1;
return y;
}
float g(float x)
{float y;
y=3*pow(x,2)-1;
return y;
}
main()
{float x0,x1,b,e,a,x;
int k;
scanf("%f,%f",&x0,&a);
k=0;
if(fabs(g(x0))<a)
{printf("error");
}else
e=1.00;
while(fabs(f(x0))>=a)
{do{x1=x0-e*f(x0)/g(x0);
b=x0;
x0=x1;
e=e/2.00;
}while(fabs(f(x1))>=fabs(f(b)));
k+=1;e=1.00;
}
x=x0;
printf("%d,%f",k,x);
}
（17）秦九韶
#include<stdio.h>
#include<math.h>
#define N 100
main()
{int m,i,j;
double sum,x,a[N];
printf("input the x and the top exp:");
scanf("%lf,%d",&x,&m);
printf("/ninput the a[i](0<i<100):");
for(i=0;i<=m;i++)
{scanf("%lf",&t);
a[i]=t;
}
sum=a[m];
for(j=m-1;j>=0;j--)
sum=x*sum+a[j];
printf("/nthe result is:%lf/n",sum);
}
（18）三对角线追赶法
#include<stdio.h>
#include<math.h>
main()
{int i,j,k,n;
float d[10][10],g[10],a[10],b[10],c[10],x[10],y[10],f[10];
printf("the top exp is ");
scanf("%d",&n);
scanf("%f,%f,%f,%f",&d[0][0],&d[0][1],&d[n-1][n-2],&d[n-1][n-1]);
for(i=1;i<n-1;i++)
for(j=i-1;j<=i+1;j++)
scanf("%f",&d[i][j]);
for(i=0;i<n;i++)
scanf("%f",&g[i]);
for(i=1;i<n-1;i++)
a[i]=d[i][i-1];
for(i=0;i<n;i++)
b[i]=d[i][i];
for(i=0;i<n-1;i++)
c[i]=d[i][i+1];
f[0]=c[0]/b[0];
for(k=1;k<n-1;k++)
f[k]=c[k]/(b[k]-a[k]*f[k-1]);
y[0]=g[0]/b[0];
for(i=1;i<n;i++)
y[i]=(g[i]-a[i]*y[i-1])/(b[i]-a[i]*f[i-1]);
x[n-1]=y[n-1];
for(i=n-2;i>=0;i--)
x[i]=y[i]-f[i]*x[i+1];
for(i=0;i<n;i++)
printf("%f/n",x[i]);
}
（19）系统
#include <stdio.h>
int main(void)
{
int a,b,CHECK;
printf("Input two nonzero intergers:");
scanf("%d,%d",&a,&b);
printf("%s%4d/n%s%4d/n%s%4d/n%s%4d/n%s%4d/n",
" a=",a,
" b=",b,
" a/b=",a/b,
" a%b=",a%b,
" CHECK=",(a/b)*b+a%b-a
);
return 0;
}
（20）弦割法
#include<stdio.h>
#include<math.h>
float f(float x)
{float y;
y=pow(x,3)-x-1;
return y;
}
main()
{float x0,x1,x2,a;
int k,N;
k=0;
scanf("%f,%f,%f,%d",&x0,&x1,&a,&N);
for(k=0;k<N;k++)
{x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0));
if(fabs(x2-x1)<a)
{printf("%d,%f",k,x2);
break;
}
else x0=x1;
x1=x2;
}
if(k>=N)
{printf("error");
}
}
（21）雅克比迭代
#include<stdio.h>
#include<math.h>
main()
{float a[10][10],b[10],x[10],y[10],e,sum,c;
int i,j,n,l;
printf("The top exp is ");
scanf("%d",&n);
printf("Now input array A/n");
for(i=0;i<n;i++)
for(j=0;j<n;j++)
scanf("%f",&a[i][j]);
printf("Now input array B first and then array X/n");
for(i=0;i<n;i++)
scanf("%f,%f",&b[i],&x[i]);
printf("Now input e/n");
scanf("%f",&e);
l=0;
do{for(i=0;i<n;i++)
{sum=0;
for(j=0;j<i;j++)
sum+=a[i][j]*x[j];
for(j=i+1;j<n;j++)
sum+=a[i][j]*x[j];
y[i]=(b[i]-sum)/a[i][i];
l+=1;
}
c=fabs(x[0]-y[0]);
for(i=0;i<n;i++)
if(c<fabs(x[i]-y[i]))
c=fabs(x[i]-y[i]);
for(i=0;i<n;i++)
x[i]=y[i];
}while(c<e);
printf("%d/n",l);
for(i=0;i<n;i++)
printf("%f/n",y[i]);
}
（22）综合测评
#include<stdio.h>
#include<math.h>
float g(float a,float b,float c,float d,float e,float f)
{float y;
y=(a*36.0+b*35.0+c*65.0+d*48.0+e*21.0+f*18.0)/223.0;
return y;
}
main()
{float a,b,c,d,e,f,num[5],max;
int i,j,h,n;
for(i=1;i<=4;i++)
{printf("input your marks:/n");
scanf("%f%f%f%f%f%f",&a,&b,&c,&d,&e,&f);
printf("your average mark:");
printf("[%d]/n",i);
printf("%6.3f/n",g(a,b,c,d,e,f));
num[i]=g(a,b,c,d,e,f);
}
for(j=1;j<=4;j++)
{max=num[1];
h=1;
for(i=1;i<=4;i++)
if(num[i]>=max)
{max=num[i];
h=i;
}
printf("%6.4f,%d/n",max,h);
num[h]=0;
}
}
2、运筹学
（1）BRANCH：分枝定界算法
#include <stdio.h>
#define len sizeof(struct node)
typedef struct node{
float bound;
int staus[50];
struct node *next;
}node;
int item[50],wl,n,state[50];
float value[50],weight[50],max_value,ratio[50];
void dele(node *father,node *current){
if(current->next==NULL)
{father->next=NULL;
return;
}
father->next=current->next;
}
void init(node *father,node *son){
int i;
father->next=son;
for(i=0;i<n;i++)
son->staus[i]=0;
son->next=NULL;
}
void branch(){
int i,t,j;
float diff,sum=0,sum_value=0;
node *head,*sonbrother,*father,*son,*prenode,*p,*q;
head=prenode=(node *)malloc(len);
father=(node *)malloc(len);
init(prenode,father);
father->bound=32768;
while(head->next!=NULL)
{
/*1*/ son=(node *)malloc(len);
init(father,son);
for(i=0;i<n&&father->staus[i]!=0;i++)
son->staus[i]=father->staus[i];
t=i;
son->staus[t]=-(t+1);
sum=0;
sum_value=0;
for(j=0;j<t+1&&son->staus[j]!=0;j++)
if(son->staus[j]>0)
{sum=sum+weight[item[j]];
sum_value=sum_value+value[item[j]];
}
while(sum!=wl&&son->staus[n-1]==0)
{diff=wl-(sum+weight[item[j]]);
if(diff>=0)
{sum=sum+weight[item[j]];
sum_value=sum_value+value[item[j]];
}
else
{sum=wl;
sum_value=sum_value+(1+diff/weight[item[j]])*value[item[j]];
}
j++;
}
son->bound=sum_value;
/*2*/ sonbrother=(node *)malloc(len);
init(son,sonbrother);
for(i=0;i<t;i++)
sonbrother->staus[i]=father->staus[i];
sonbrother->staus[t]=t+1;
sum=0;
sum_value=0;
for(j=0;j<t+1&&sonbrother->staus[j]!=0;j++)
if(sonbrother->staus[j]>0)
{sum=sum+weight[item[j]];
sum_value=sum_value+value[item[j]];
}
if(sum>wl)
{sonbrother->bound=-32768;
dele(son,sonbrother);
}
else
{while(sum!=wl&&sonbrother->staus[n-1]==0)
{diff=wl-(sum+weight[item[j]]);
if(diff>=0)
{sum=sum+weight[item[j]];
sum_value=sum_value+value[item[j]];
}
else
{sum=wl;
sum_value=sum_value+(1+diff/weight[item[j]])*value[item[j]];
}
j++;
}
sonbrother->bound=sum_value;
}
dele(prenode,father);
father=prenode->next;
if(son->staus[n-1]!=0)
{if(son->next!=NULL)
{max_value=sonbrother->bound;
for(i=0;i<n;i++)
state[i]=sonbrother->staus[i];
dele(son,sonbrother);
dele(prenode,father);
father=prenode->next;
}
else
{max_value=son->bound;
for(i=0;i<n;i++)
state[i]=son->staus[i];
dele(prenode,father);
}
q=head;
p=head->next;
while((p!=NULL)&&(p->bound<=max_value))
{dele(q,p);
p=q->next;
}
if(p!=NULL)
{prenode=q;
father=p;
}
else
return;
}
else
if(father->next!=NULL)
{prenode=prenode->next;
father=father->next;
}
}
return;
}
int getmin(){
int i;
float amin=weight[0];
for(i=1;i<n;i++)
if(amin>weight[i])
amin=weight[i];
return amin;
}
void sort(){
int i,j,exchange=1;
float temp1,temp2;
for(i=0;i<n;i++)
ratio[i]=value[i]/weight[i];
for(j=n-1;j>=0&&exchange==1;j--)
{exchange=0;
for(i=0;i<j;i++)
if(ratio[i+1]>ratio[i])
{exchange=1;
temp1=ratio[i+1];ratio[i+1]=ratio[i];ratio[i]=temp1;
temp2=item[i+1];item[i+1]=item[i];item[i]=temp2;
}
}
}
void main(){
int i,j;
float sum=0;
clrscr();
printf("/tWelcome to the BRANCH_BOUND system!/n");
printf("number of the materials=? ");
scanf("%d",&n);
printf("maximun weigh of the problem=? ");
scanf("%d",&wl);
for(i=0;i<n;i++)
{item[i]=i;
printf("/n*******************/n");
printf("input item%d data!/n",i+1);
printf("*******************/n");
printf("weight %d=? ",i+1);
scanf("%f",&weight[i]);
printf("value %d=? ",i+1);
scanf("%f",&value[i]);
}
if((getmin())>wl)
{printf("/nThere is no solution of the problem!");
exit(0);
}
for(i=0;i<n;i++)
sum=sum+weight[i];
if(sum<=wl)
{printf("/nAll the materials can be loaded!");
exit(0);
}
sort();
branch();
printf("/n/n/nThe maximum value of the materials is %f ",max_value);
printf("/nincluding the following materials/n");
sum=0;
for(i=0;i<n;i++)
if(state[i]>0)
{sum=sum+weight[item[i]];
printf("%d/n",item[i]+1);
}
printf("/nThe weight of the materials is %f ",sum);
getch();
}
（2）CHAIN：马尔可夫链算法
#include <stdio.h>
#include <math.h>
double a[10][10];
void Guass(int n){
int i,j,k;
double t;
for(k=0;k<n-1;k++)
{t=a[k][k];
for(j=k;j<n;j++)
a[k][j]=a[k][j]/t;
for(i=0;i<n-1;i++)
if(i!=k)
{t=a[i][k]/a[k][k];
for(j=k;j<n;j++)
a[i][j]=a[i][j]-a[k][j]*t;
}
}
return;
}
void chain(){
static double p[10][10],pr[10],diff,table[100][10],pnew[10][10],ptemp[10][10],temp[10],exr[10][10];
int n,i,j,k,s,m,found,inr,inc;
printf("Welcome to the MARKOV CHAIN ANALYSIS system!/n/n");
printf("how many states =? ");
scanf("%d",&n);
printf("/nthe steady transmit possibility of step 1 ?/n");
for(i=0;i<n;i++)
for(j=0;j<n;j++)
scanf("%lf",&p[i][j]);
printf("/nthe initiate state of step 1 ?/n");
for(i=0;i<n;i++)
scanf("%lf",&pr[i]);
for(i=0;i<n;i++)
for(j=0;j<n;j++)
pnew[i][j]=p[i][j];
for(i=0;i<n;i++)
table[0][i]=pr[i];
printf("/n step 1");
for(i=0;i<n;i++)
{printf("/n");
for(j=0;j<n;j++)
printf("%f/t",p[i][j]);
}
printf("/n");
for(k=2;k<100;k++)
{for(j=0;j<n;j++)
{temp[j]=0;
for(i=0;i<n;i++)
temp[j]=temp[j]+pr[i]*pnew[i][j];
}
for(i=0;i<n;i++)
table[k-1][i]=temp[i];
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{ptemp[i][j]=0;
for(m=0;m<n;m++)
ptemp[i][j]=ptemp[i][j]+p[i][m]*pnew[m][j];
}
for(i=0;i<n;i++)
for(j=0;j<n;j++)
pnew[i][j]=ptemp[i][j];
for(j=0;j<n;j++)
{for(i=0;i<n-1;i++)
{
for(m=i+1;m<n;m++)
{diff=pnew[i][j]-pnew[m][j];
if(diff<0)
diff=-diff;
if(diff>0.001)
{found=0;
break;
}
found=1;
}
if(diff>0.001) break;
}
if(diff>0.0001) break;
}
if(found==0)
{if(k%5==0)
{printf("/n step %d",k);
for(i=0;i<n;i++)
{printf("/n");
for(j=0;j<n;j++)
printf("%f/t",pnew[i][j]);
}
getch();
}
if(k>=100)
{printf("/nsteady_state probability have not been detained in 100");
return;
}
}
else
{printf("/nstep %d",k);
for(i=0;i<n;i++)
{printf("/n");
for(j=0;j<n;j++)
printf("%f/t",pnew[i][j]);
}
break;
}
}
for(j=0;j<n;j++)
{temp[j]=0;
for(i=0;i<n;i++)
temp[j]=temp[j]+pr[i]*pnew[i][j];
}
for(i=0;i<n;i++)
table[k][i]=temp[i];
printf("/nThe steady-state probability of being in /n");
for(j=0;j<n;j++)
printf("state %d is %f/n",j,pnew[n-1][j]);
printf("probability of being in state/n/n");
for(i=0;i<=k;i++)
{printf("%d",i);
for(j=0;j<n;j++)
printf("/t%f",table[i][j]);
printf("/n");
if(i%10==0) getch();
}
for(s=0;s<n;s++)
{ inr=0;
for(j=0;j<n;j++)
if(j==s)
continue;
else
{inc=0;
for(i=0;i<n;i++)
if(i==s)
continue;
else
{
a[inr][inc]=-p[j][i];
if(j==i)
a[inr][inc]=1+a[inr][inc];
inc++;
}
inr++;
}
for(i=0;i<n-1;i++)
a[i][n-1]=1;
Guass(n);
i=0;
for(j=0;j<n;j++)
if(j!=s)
exr[j][s]=a[i++][n-1];
else
exr[j][s]=1/pnew[n-1][s];
}
printf("/nTable of expected first passage times and recurrence times");
for(i=0;i<n;i++)
{printf("/n%d",i);
for(j=0;j<n;j++)
printf("/t%f",exr[i][j]);
}
}
void main(){
clrscr();
chain();
getch();
}
（3）DECISION：贝叶斯决策方法
#include <stdio.h>
#include <math.h>
#define pi 3.14159
#define p(x,t) exp(-(x-t)*(x-t)/20)/sqrt(20*pi)
void decision(){
int i,j,type,m,n,flag,state[5],index;
float xx,a[5][5],p[5],e[5],sum,decision;
printf("Welcome to the DECISION_STSTEM!");
printf("/n/ntype of the problem,max(key '0')or min(key '1')? ");
scanf("%d",&type);
printf("type of the decision,without data(key'0')or with data(key'1')? ");
scanf("%d",&flag);
printf("number of the actions ");
scanf("%d",&m);
printf("number of the nature states ");
scanf("%d",&n);
for(j=0;j<n;j++)
printf("/tt%d",j+1);
printf("/n");
for(i=0;i<m;i++)
{printf("A%d/t",i+1);
for(j=0;j<n;j++)
scanf("%f",&a[i][j]);
}
printf("probability of the natural states=?/n");
for(j=0;j<n;j++)
scanf("%f",&p[j]);
if(flag)
{printf("given data=? ");
scanf("%f",&xx);
printf("states of nature=?/n");
for(j=0;j<n;j++)
scanf("%d",&state[j]);
sum=0;
for(j=0;j<n;j++)
sum+=p[j]*p(xx,state[j]);
for(j=0;j<n;j++)
p[j]=p[j]*p(xx,state[j])/sum;
}
for(i=0;i<m;i++)
{sum=0;
for(j=0;j<n;j++)
sum+=p[j]*a[i][j];
e[i]=sum;
}
decision=e[0];
index=1;
if(type==0)
for(i=1;i<m;i++)
{ if(decision<e[i])
{decision=e[i];
index=i+1;
}
}
else
for(i=0;i<m;i++)
if(decision>e[i])
{decision=e[i];
index=i+1;
}
printf("/n**********/n");
printf("Results:");
printf("/n**********/n");
printf("expected loss for each course of action based on prior distribution/n");
for(i=0;i<m;i++)
printf("%f/t",e[i]);
printf("/n/nThe optimum expected loss is %f ",decision);
printf("/n/nchoose action A( %d )",index);
return;
}
void main(){
clrscr();
decision();
getch();
}
（4）dp_invest：动态规划的投资问题
#include <stdio.h>
int istar[10];
float dp_invest(int N,int K){
int i,j,sum,z,d[10][50];
float g[10][50],f[10][50];
printf("/nThe return function values as follows!/n");
for(j=0;j<K+1;j++)
printf("/t%d",j);
for(i=0;i<N;i++)
{printf("/ntask%d:/t",i+1);
for(j=0;j<K+1;j++)
scanf("%f",&g[i][j]);
}
for(i=0;i<N;i++)
for(j=0;j<K+1;j++)
{f[i][j]=0;d[i][j]=0;}
for(j=0;j<K+1;j++)
{f[N-1][j]=g[N-1][j];
d[N-1][j]=j;
}
for(i=N-2;i>=0;i--)
for(j=1;j<K+1;j++)
{f[i][j]=g[i][0]+f[i+1][j];
d[i][j]=0;
for(z=1;z<=j;z++)
if((g[i][z]+f[i+1][j-z])>f[i][j])
{f[i][j]=g[i][z]+f[i+1][j-z];
d[i][j]=z;
}
}
istar[0]=d[0][K];
for(i=1;i<N;i++)
{sum=0;
for(j=0;j<=i-1;j++)
sum=sum+istar[j];
istar[i]=d[i][K-sum];
}
return f[0][K];
}
void main(){
int i,N,K;
clrscr();
printf("WELCOME TO THE DYNAMIC_INVEST SYSTEM!/n/n");
printf("How many tasks ? ");
scanf("%d",&N);
printf("/nHow many units of materials ? ");
scanf("%d",&K);
printf("/nThe optimal return is %f",dp_invest(N,K));
for(i=0;i<N;i++)
printf("/nThe optimal amout to invest in task %d is %d",i+1,istar[i]);
getch();
}
（5）dp_plan：生产计划算法
#include <stdio.h>
#define pc(j) 20+5*j
#define e(j) j
void dp_plan(){
int i,j,k,sum,limit,n,io,max_stock,max_production,d[10],x[10][100],z,z_maxlimit,z_minlimit,xstar[10];
float f[10][100],temp;
printf("WELCOME TO THE DYNIMIC SYSTEM!/n/nperiods of the inventory =? ");
scanf("%d",&n);
printf("/nthe maximum stocks of each period=? ");
scanf("%d",&max_stock);
printf("/nthe maximum production of each period=? ");
scanf("%d",&max_production);
printf("/nThe stocks of the first period=? ");
scanf("%d",&io);
for(i=0;i<n;i++)
{printf("Demand for period %d is ",i+1);
scanf("%d",&d[i]);
}
for(k=0;k<d[n-1];k++)
{x[n-1][k]=d[n-1]-k;
f[n-1][k]=pc(x[n-1][k]);
}
f[n-1][d[n-1]]=0;
x[n-1][d[n-1]]=0;
for(i=n-2;i>=0;i--)
{sum=0;
for(j=i;j<n;j++)
sum+=d[j];
if(sum<max_stock)
limit=sum;
else
limit=max_stock;
for(k=0;k<=limit;k++)
{if(d[i]-k>0)
{z_minlimit=d[i]-k;
f[i][k]=pc(z_minlimit)+e(0)+f[i+1][0];
x[i][k]=z_minlimit;
}
else
{z_minlimit=0;
f[i][k]=e(k-d[i])+f[i+1][k-d[i]];
x[i][k]=0;
}
if(sum-k>max_stock+d[i]-k)
if(max_stock+d[i]-k>max_production)
z_maxlimit=max_production;
else
z_maxlimit=max_stock+d[i]-k;
else
if(sum-k>max_production)
z_maxlimit=max_production;
else
z_maxlimit=sum-k;
for(z=z_minlimit;z<=z_maxlimit;z++)
{temp=pc(z)+e(k+z-d[i])+f[i+1][k+z-d[i]];
if(f[i][k]>temp)
{f[i][k]=temp;
x[i][k]=z;
}
}
}
}
/* for(i=0;i<n;i++)
{printf("/n/nthe period %d/n",i+1);
for(j=0;j<=5;j++)
printf("/t%f--->%d/t",f[i][j],x[i][j]);
getch();
}*/
printf("/n/nThe minimum policy cost for the %d periods is %f",n,f[0][io]);
xstar[0]=x[0][io];
j=io;
printf("/n/nThe optimal amount to produce in period 1 is %d",xstar[0]);
for(i=1;i<n;i++)
{xstar[i]=x[i][xstar[i-1]-d[i-1]+j];
printf("/nThe optimal amount to produce in period %d is %d",i+1,xstar[i]);
j=xstar[i-1]-d[i-1]+j;
}
}
void main(){
clrscr();
dp_plan();
getch();
}
（6）linear：线性规划单纯形算法
#include <stdio.h>
int K,M,N,Q=100,Type,Get,Let,Et,Code[50],XB[50],IA,IAA[50],Indexg,Indexl,Indexe;
float Sum,A[50][50],B[50],C[50];
void initiate();
void solve();
void main(){
int i,j;
clrscr();
/****** input data ******/
printf("/nWelcome to the linear programming solution!/n/n********/nNotice 1!/n********/nIf the type of objective equation is 'max' ,please enter '1'!/nElse please enter '0'!/n/n/n");
printf("********/nNotice 2!/n********/nDefine the type of subjective equation as following:/n'<='is equal to '0'!/n'>=' is equal to '1'!/n'='is equal to '2'!/n");
printf("/n/nPlease input the coefficients or the constants:/n");
printf("THe number of subjective equations ? ");
scanf("%d",&M);
printf("THe number of variables ? ");
scanf("%d",&K);
printf("THe number of ' <= 'subjective equations ? ");
scanf("%d",&Let);
printf("THe number of ' >= 'subjective equations ? ");
scanf("%d",&Get);
printf("THe number of ' = 'subjective equations ? ");
scanf("%d",&Et);
printf("THe type of objective equation ? ");
scanf("%d",&Type);
N=K+Let+Et+2*Get;
for(i=0;i<M;i++)
{printf("Please input %d's equation:/n",i+1);
printf("type ? ");
scanf("%d",&Code[i]);
printf("constant ? ");
scanf("%f",&B[i]);
for(j=0;j<K;j++)
{printf("coefficient ? ");
scanf("%f",&A[i][j]);
}
}
printf("Plese input the constants of the object equation:/n");
for(j=0;j<K;j++)
scanf("%f",&C[j]);
printf("/n/n*************************************/nPlease check the data you just input!/n************************************* /n");
getch();
if(Type)
printf("The type of the object equation is 'max'/n");
else
printf("The type of the object equation is 'min'/n");
printf("/nThe object equation is :/n");
for(j=0;j<K;j++)
{printf("(%f)X%d ",C[j],j+1);
if(j!=K-1)
printf("+");
}
printf("/n/nThe suject equation is :");
for(i=0;i<M;i++)
{printf("/nNumber %d suject equation is: %d %f ",i+1,Code[i],B[i]);
for(j=0;j<K;j++)
{printf("(%f)X%d ",A[i][j],j+1);
if(j!=K-1)
printf("+");
}
}
/****** initiate data ******/
initiate();
solve();
if(!Type)
A[M][N]=-A[M][N];
printf("/n/nThe optimal value of the original objective function is: %f ",A[M][N]);
getch();
}
/****** initiate variables function ******/
void initiate(){
int i,j;
Indexg=K;
Indexl=Indexg+Get;
Indexe=Indexl+Let;;
for(i=0;i<M+1;i++)
for(j=K;j<N+1;j++)
A[i][j]=0;
for(i=0;i<M;i++)
A[i][N]=B[i];
for(i=0;i<M;i++)
switch(Code[i])
{ case 0: {XB[i]=Indexl;
A[i][Indexl++]=1;
break;
};
case 1: {XB[i]=Indexe;
IAA[IA++]=i;
A[i][Indexe++]=1;
A[i][Indexg++]=-1;
break;
};
case 2: {XB[i]=Indexe;
IAA[IA++]=i;
A[i][Indexe++]=1;
break;
};
}
for(j=0;j<K;j++)
if(Type)
A[M][j]=-C[j];
else
A[M][j]=C[j];
for(j=K;j<=N;j++)
A[M][j]=0;
for(j=K+Get+Let;j<N;j++)
A[M][j]=Q;
Sum=0;
for(j=0;j<=N;j++)
{Sum=0;
for(i=0;i<IA;i++)
Sum=Sum+A[IAA[i]][j];
A[M][j]=A[M][j]-Sum*Q;
}
return;
}
/****** process data function ******/
void solve(){
int i,j,mark=1,minus,minusmark,basic=0,divide,dividemark;
float H,P;
while(1)
{mark=0;
minusmark=0;
minus=0;divide=1000;
dividemark=0;
printf("/n/nBasic solution %d is/n",++basic);
for(i=0;i<M;i++)
printf("Basic variable %d = X( %d )= %f/n",i+1,XB[i]+1,A[i][N]);
printf("/nCurrent value of the object equation is: %f/n",A[M][N]);
getch();
for(j=0;j<N;j++)
{
if(A[M][j]<-6e-8)
mark++;
if(A[M][j]<minus)
{minus=A[M][j];
minusmark=j;
}
}
if(mark==0)
break;
for(i=0;i<M;i++)
{if(A[i][minusmark]==0)
continue;
if(A[i][N]/A[i][minusmark]<=0)
continue;
if(A[i][N]/A[i][minusmark]<divide)
{divide=A[i][N]/A[i][minusmark];
dividemark=i;
}
}
XB[dividemark]=minusmark;
if(divide<0)
printf("There is no solution because of no boundary!");
P=A[dividemark][minusmark];
for(j=0;j<N+1;j++)
A[dividemark][j]=A[dividemark][j]/P;
for(i=0;i<M+1;i++)
{H=A[i][minusmark];
if(i==dividemark)
continue;
for(j=0;j<N+1;j++)
A[i][j]=A[i][j]-H*A[dividemark][j];
}
}
printf("/n/n*************************************/nThe last basic solution is optimal!/n************************************* /n");
return;
}
扫一扫
2388
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
