[原创]单纯形法完全c语言程序,能运行

sprine

单纯形法完全c语言程序，每步都有说明，不足之处希望指出。

1. #include "math.h"
   
2. #include "stdio.h"
   
3. #define N 2
   
4. 
   
5. void paixu(p,n)
   
6. int n;
   
7. double p[];
   
8. { int m,k,j,i;
   
9. double d;
   
10. k=0; m=n-1;
    
11. while (k<m)
    
12. { j=m-1; m=0;
    
13. for (i=k; i<=j; i++)
    
14. if (p>p[i+1])
    
15. { d=p; p=p[i+1]; p[i+1]=d; m=i;}
    
16. j=k+1; k=0;
    
17. for (i=m; i>=j; i--)
    
18. if (p[i-1]>p)
    
19. { d=p; p=p[i-1]; p[i-1]=d; k=i;}
    
20. }
    
21. return;
    
22. }
    
23. double mubiao(double *x)
    
24. { double y;
    
25. y=x[1]-x[0]*x[0]; y=100.0*y*y;
    
26. y=y+(1.0-x[0])*(1.0-x[0]);
    
27. return(y);
    
28. }
    
29. 
    
30. main()
    
31. { int i,j,k,l,m=0;
    
32. double c,xx[N+1][N],f0[N+1],f[N+1],x0[N]={1.2,1},x1[N],s=0.0;
    
33. double a,b;
    
34. double xa[N],xb[N],xc[N],xe[N],xw[N],xr[N],xo[N];
    
35. double fr,fe,fw,fc,fo;
    
36. double aef=1.0,r=1.0,eps1=1.0e-30,eps2=1.0e-30,bt=0.5,rou=0.5;
    
37. c=1.0;
    
38. b=(c/(N*sqrt(2)))*(sqrt(N+1)-1);
    
39. a=b+c/sqrt(2);
    
40. //printf("a=%13.7e b=%13.7e ",a,b);
    
41. //printf("\n");
    
42. //给xx[N][N+1]赋值,每一行构成单纯形的一个定点
    
43. //***********************
    
44. for(i=0;i<N;i++)
    
45. xx[0]=x0;
    
46. for(i=1;i<N+1;i++)
    
47. for(j=0;j<N;j++)
    
48. {if(j==i-1)
    
49. xx[j]=x0[j]+a;
    
50. else
    
51. xx[j]=x0[j]+b;
    
52. }
    
53. for (i=0;i<N+1;i++)
    
54. {for (j=0;j<N;j++)
    
55. printf("xx[%d][%d]%13.7e ",i,j,xx[j]);
    
56. printf("\n");
    
57. }
    
58. loop1:
    
59. //求单纯形的每个定点的函数值f0,f和x1是过渡数组
    
60. printf("\n");
    
61. printf("\n");
    
62. for(i=0;i<N+1;i++)
    
63. {for(j=0;j<N;j++)
    
64. x1[j]=xx[j];
    
65. f0=mubiao(x1);
    
66. f=mubiao(x1);
    
67. printf("f0[%d]=%13.7e f[%d]=%13.7e\n",i,f0,i,f);
    
68. }
    
69. printf("\n");
    
70. //比较f的大小，f[0]是最小值，f[N]是最大值
    
71. paixu(f,N+1);
    
72. for(i=N;i>=0;i--)
    
73. printf("f[%d]=%13.7e \n",i,f);
    
74. //找最好点和最坏点分别是哪一个点，即xx[][]的行数
    
75. for(i=0;i<N+1;i++)
    
76. {if(f0==f[0])
    
77. k=i;
    
78. if(f0==f[N])
    
79. l=i;
    
80. }
    
81. printf("最好点k=%d\n",k);
    
82. printf("最坏点l=%d\n",l);
    
83. //终止判断条件
    
84. printf("f[N]-f[0]=%13.7e \n",f[N]-f[0]);
    
85. if((f[N]-f[0])<eps1+eps2*fabs(f[N]))
    
86. {printf("迭代次数m=%d\n",m);
    
87. for(j=0;j<N;j++)
    
88. printf("optx[%d]=%13.7e\n",j,xx[k][j]);
    
89. printf("fmin=%13.7e\n",f[0]);
    
90. }
    
91. else
    
92. {
    
93. m=m+1;
    
94. //把xx[][]中最好点移到第一行，最坏点移到最后一行
    
95. for(j=0;j<N;j++)
    
96. {xb[j]=xx[k][j];
    
97. xx[k][j]=xx[0][j];
    
98. xx[0][j]=xb[j];
    
99. //
    
100. xw[j]=xx[l][j];
     
101. xx[l][j]=xx[N][j];
     
102. xx[N][j]=xw[j];
     
103. }
     
104. for (i=0;i<N+1;i++)
     
105. {for (j=0;j<N;j++)
     
106. printf("xx[%d][%d]=%13.7e ",i,j,xx[j]);
     
107. printf("\n");
     
108. }
     
109. //求除最坏点f[N]外其余点的中点xc[]
     
110. for(i=0;i<N;i++)
     
111. xa=0;
     
112. for(j=0;j<N;j++)
     
113. {{for(i=0;i<N;i++)
     
114. 
     
115. xa[j]=xa[j]+xx[j];}
     
116. xa[j]=xa[j]/N;
     
117. }
     
118. for(i=0;i<N;i++)
     
119. printf("xa[%d]=%13.7e xb[%d]=%13.7e xw[%d]=%13.7e \n",i,xa,i,xb,i,xw);
     
120. //求xw[N]的反射点xr[N];
     
121. for(i=0;i<N;i++)
     
122. {xr=xa+aef*(xa-xw);
     
123. printf("xr[%d]=%13.7e ",i,xr);
     
124. }
     
125. printf("\n");
     
126. //求xr[N]的函数值fr
     
127. fr=mubiao(xr);
     
128. printf("fr=%13.7e \n",fr);
     
129. //判断xr与xb的好坏
     
130. if(fr<=f[0])
     
131. {for(i=0;i<N;i++)
     
132. {xe=xr+r*(xr-xa);
     
133. //printf("xe[%d]=%13.7e ",i,xe);
     
134. }
     
135. printf("\n");
     
136. fe=mubiao(xe);
     
137. if(fe<=f[0])
     
138. for(j=0;j<N;j++)
     
139. xx[N][j]=xe[j];
     
140. else
     
141. for(j=0;j<N;j++)
     
142. xx[N][j]=xr[j];
     
143. goto loop1;
     
144. }
     
145. else
     
146. {
     
147. fw=f[N];
     
148. if(fr>=fw)
     
149. {for(i=0;i<N;i++)
     
150. xc=xa-bt*(xa-xw);
     
151. fc=mubiao(xc);
     
152. if(fc>=fw)
     
153. {for(i=1;i<N+1;i++)
     
154. for(j=0;j<N;j++)
     
155. xx[j]=xx[0][j]-rou*(xx[j]-xx[0][j]);
     
156. goto loop1;
     
157. }
     
158. else
     
159. {for(j=0;j<N;j++)
     
160. xx[N][j]=xc[j];
     
161. goto loop1;
     
162. }
     
163. }
     
164. else
     
165. {if(fr>=fe)
     
166. {
     
167. for(i=0;i<N;i++)
     
168. xo=xa+bt*(xa-xw);
     
169. fo=mubiao(xo);
     
170. if(fo>=fr)
     
171. {for(i=1;i<N+1;i++)
     
172. for(j=0;j<N;j++)
     
173. xx[j]=xx[0][j]-rou*(xx[j]-xx[0][j]);
     
174. goto loop1;
     
175. }
     
176. else
     
177. {for(j=0;j<N;j++)
     
178. xx[N][j]=xo[j];
     
179. goto loop1;
     
180. }
     
181. }
     
182. else
     
183. {for(j=0;j<N;j++)
     
184. xx[N][j]=xr[j];
     
185. goto loop1;
     
186. }
     
187. }
     
188. }
     
189. }
     
190. }

复制代码

[ 本帖最后由 风花雪月 于 2006-10-18 08:23 编辑 ]

miaoyu0807

你的程序中的各个变量都代表什么意思？
谢谢！！！

zhuwy16

好像不能够运行啊

sprine

原帖由 zhuwy16 于 2006-7-5 18:10 发表
好像不能够运行啊

我是在VC6.0中编译通过的，再试试吧！！

winter33123152

我以前也编过一个,不过收敛性不怎么好.有时间试试楼主的怎么样.