黄金分割法-机械优化设计程序
- 格式:doc
- 大小:41.00 KB
- 文档页数:6
一、实验目的本次实验旨在通过计算机编程,加深对机械优化设计方法的理解,掌握常用的优化算法,并能够利用计算机解决实际问题。
二、实验内容1. 黄金分割法(1)实验原理黄金分割法是一种常用的优化算法,适用于一元函数的极值求解。
其基本原理是:在给定初始区间内,通过迭代计算,逐步缩小搜索区间,直到满足收敛条件。
(2)实验步骤① 设计实验程序,实现黄金分割法的基本算法。
② 编写函数,用于计算一元函数的值。
③ 设置初始区间和收敛精度。
④ 迭代计算,更新搜索区间。
⑤ 判断是否满足收敛条件,若满足则输出结果,否则继续迭代。
(3)实验结果通过编程实现黄金分割法,求解函数f(x) = x^3 - 6x^2 + 9x + 1在区间[0, 10]内的极小值。
实验结果显示,该函数在区间[0, 10]内的极小值为1,且收敛精度达到0.001。
2. 牛顿法(1)实验原理牛顿法是一种求解非线性方程组的优化算法,其基本原理是:利用函数的导数信息,逐步逼近函数的极值点。
(2)实验步骤① 设计实验程序,实现牛顿法的基本算法。
② 编写函数,用于计算一元函数及其导数。
③ 设置初始值和收敛精度。
④ 迭代计算,更新函数的近似值。
⑤ 判断是否满足收敛条件,若满足则输出结果,否则继续迭代。
(3)实验结果通过编程实现牛顿法,求解函数f(x) = x^3 - 6x^2 + 9x + 1在区间[0, 10]内的极小值。
实验结果显示,该函数在区间[0, 10]内的极小值为1,且收敛精度达到0.001。
3. 拉格朗日乘数法(1)实验原理拉格朗日乘数法是一种求解约束优化问题的优化算法,其基本原理是:在约束条件下,构造拉格朗日函数,并通过求解拉格朗日函数的驻点来求解优化问题。
(2)实验步骤① 设计实验程序,实现拉格朗日乘数法的基本算法。
② 编写函数,用于计算目标函数、约束函数及其导数。
③ 设置初始值和收敛精度。
④ 迭代计算,更新拉格朗日乘数和约束变量的近似值。
黄金分割法机械优化设计在现代工程设计领域,机械优化设计是一项非常重要的任务。
通过对机械系统进行分析和优化,可以提高其性能和效率,节约资源并延长使用寿命。
黄金分割法是一种常用的优化设计方法,它基于黄金分割比的原理,通过寻找最佳设计参数来改进机械系统的性能。
本文将介绍黄金分割法机械优化设计的原理、方法和应用。
一、黄金分割法的原理黄金分割法源自于数学中的黄金分割比,即0.618,也称为费波那契数。
它是指将一条线段分割为两部分,使较长部分与整体的长度之比等于较短部分与较长部分之比。
黄金分割法的原理是将这一比例应用于机械设计中,以找到最佳的设计参数。
二、黄金分割法机械优化设计的方法1. 确定优化目标:在机械优化设计中,首先需要明确具体的优化目标。
比如,改善机械系统的运行效率、减少能源消耗或提高产品质量等。
2. 确定设计参数:根据机械系统的特性和优化目标,确定需要进行优化的设计参数。
这些参数可以是机械结构的尺寸、材料的选择或运行参数等。
3. 建立优化模型:根据设计参数,建立机械系统的优化模型。
模型可以是数学模型、仿真模型或实验模型,根据具体情况选择。
4. 寻找最佳设计参数:利用黄金分割法进行参数优化。
通过分割设计参数范围,并根据黄金分割比的原理,逐步缩小搜索范围,最终找到最佳设计参数。
5. 评估和验证:对优化得到的设计参数进行评估和验证。
可以通过数值模拟、物理实验或现场测试等方法,验证优化结果是否满足设计要求。
三、黄金分割法机械优化设计的应用黄金分割法机械优化设计在各行业都有广泛的应用。
以下为几个常见的应用领域:1. 机械结构设计:对于机械结构的设计优化,黄金分割法可以帮助确定最佳的尺寸比例,提高结构的刚性和稳定性。
2. 流体力学设计:在流体力学设计中,黄金分割法可以通过优化设计参数,改善流体的流动性能,提高流体的传输效率和混合效果。
3. 电子电路设计:黄金分割法可以应用于电子电路设计中,通过优化电路元件的参数和布局来提高电路的性能和稳定性。
机械优化黄金分割法程序设计机械优化黄金分割法程序设计1.黄金分割法介绍黄金分割法适用于[a,b]区间上的任何单股函数求极小值问题,对函数除要求“单谷”外不做其他要求,甚至可以不连续。
因此,这种方法的适应面非常广。
黄金分割法也是建立在区间消去法原理基础上的试探方法,即在搜索区间[a,b]内适当插入两点a1,a2,并计算其函数值。
a1,a2将区间分成三段,应用函数的单谷性质,通过函数值大小的比较,删去其中一段,是搜索区间得以缩小。
然后再在保留下来的区间上作同样的处理,如此迭代下去,是搜索区间无限缩小,从而得到极小点的数值近似解。
1.1黄金分割法原理一维搜索是解函数极小值的方法之一,其解法思想为沿某一已知方向求目标函数的极小值点。
一维搜索的解法很多,这里主要采用黄金分割法(0.618法)。
该方法用不变的区间缩短率0.618代替斐波那契法每次不同的缩短率,从而可以看成是斐波那契法的近似,实现起来比较容易,也易于人们所接受。
黄金分割法是用于一元函数f(x)在给定初始区间[a,b]内搜索极小点α*的一种方法。
它是优化计算中的经典算法,以算法简单、收敛速度均匀、效果较好而著称,是许多优化算法的基础,但它只适用于一维区间上的凸函数[6],即只在单峰区间内才能进行一维寻优,其收敛效率较低。
其基本原理是:依照“去劣存优”原则、对称原则、以及等比收缩原则来逐步缩小搜索区间[7]。
具体步骤是:在区间[a,b]内取点:a1 ,a2 把[a,b]分为三段。
如果f(a1)>f(a2),令a=a1,a1=a2,a2=a+r*(b-a);如果f(a1)黄金分割法:黄金分割法适用于确定区间上的任何单谷函数求极小值的问题。
对函数除要求“单谷”之外没有任何其他要求。
1. 给出初始搜索区间,及收敛精度将其赋以0.618。
2. 计算a1和a2,并计算起对应的函数值y1,y2。
3. 根据期间消去法原理缩短搜索区间,为了能用原来的坐标点计算公式,需进行区间名称的代换,并在保留区间中计算一个新的试验点及其函数值。
课程设计(实验)材料(2)机械优化设计课程设计(实验)报告专业班级:设计题目: 黄金分割法程序设计学生姓名:学生学号:任课教师:2013年月日一、设计要求:基于一维搜索的试探方法思想,运用黄金分割法编写C语言程序,得到极小点的数值近似解。
已知条件:1、目标函数:f(x)=x*x+2*x2、给定搜索区间:(-3,5)二、方法原理黄金分割法适用于【a,b】区间上的任何单股函数求极小值问题,对函数除要求“单谷”外不做其他要求,甚至可以不连续 .因此 ,这种方法的适应面非常广.黄金分割法也是建立在区间消去法原理基础上的试探方法,即在搜索区间【a,b】内适当插入两点a1,a2,并计算其函数值.a1,a2将区间分成三段,应用函数的单谷性质,通过函数值大小的比较,删去其中一段,是搜索区间得以缩小。
然后再在保留下来的区间上作同样的处理,如此迭代下去,是搜索区间无限缩小,从而得到极小点的数值近似解。
三、程序清单:#include<stdio.h〉#include〈math.h〉double hanshu (double x);void main(){int k;double a,a1,a2,b,y1,y2,c,e,i,j;e=0.618,k=0;printf(”a=");scanf(”%lf",&a);printf("b=");scanf("%lf”,&b);printf(”c=");scanf(”%lf”,&c);a1=b-e*(b-a);y1=hanshu(a1);a2=a+e*(b—a);y2=hanshu(a2);printf("输出次数=%d\n a=%lf,a1=%lf, a2=%lf,b=%lf,y1=%lf, y2=%lf\n”,k,a, a1,a2,b,y1,y2);i=(b-a)/b;j=(y2—y2)/y2;while(fabs(i)〉=c||fabs(j>=c)){k++;if(y1〉=y2){a=a1,a1=a2,y1=y2,a2=a+e*(b-a),y2=hanshu(a2);}else{b=a2,a2=a1,y2=y1,a1=b-e*(b-a),y1=hanshu(a1);}printf(”输出次数=%d\n a=%lf,a1=%lf, a2=%lf, b=%lf,y1=%lf, y2=%lf\n”,k,a,a1,a2,b,y1,y2);i=(b-a)/b;j=(y2—y1)/y2;}printf("k=%lf\n”,0.5*(a+b));}double hanshu (double x){double m;m = x*x+2*x;return m;}实验结果(要求附上程序运行结果截图)五、手算过程如以下表格:。
目录一、黄金分割法二、二次插值法三、最速下降法(阶梯法)四、变尺度法五、鲍威尔法一、黄金分割法#include<stdio.h>#include<math。
h>#define r 0。
618#define f(x)x*x+2*xgolden(float,float,float);main(){float a,b,e;printf(”\n请输入区间和收敛精度:a,b,e\n”);scanf("%f,%f,%f”,&a,&b,&e);golden(a,b,e);}golden(float a,float b,float e){float y1,y2,a1,a2,A,Y;int n=0;a1=b—r*(b-a);a2=a+r*(b-a);y1=f(a1);y2=f(a2);printf(”黄金分割法的搜索过程:");do{printf("\n %d a=%f,b=%f,a1=%f,a2=%f,y1=%f,y2=%f",n,a,b,a1,a2,y1,y2);if(y1〉=y2){a=a1;a1=a2;y1=y2;a2=a+r*(b—a);y2=f(a2);}else{b=a2;a2=a1;y2=y1;a1=b-r*(b—a);y1=f(a1);}n++;}while(fabs((b—a)/b)〉=e||fabs((y2-y1)/y2)>=e);A=(a+b)/2;Y=f(A);printf("\n %d a=%f,b=%f,a1=%f,a2=%f”,n,a,b,a1,a2);printf("\n结果:\n极值点及其函数值:A=%f,Y=%f\n”, A,Y);}二、二次插值法#include"stdio.h"#include”math.h”#include”conio.h"void main(){float*area(float a1,float p,float a[3]);float f(float x);float ar,fr;float a1=10,p=0.01,e=0。
#include"stdio.h"#include"stdlib.h"#include"math.h"double objf(double x[]){double ff;ff=x[0]*x[0]+x[1]*x[1]-x[0]*x[1]-10*x[0]-4*x[1]+60;return(ff);}double gold(double a[],double b[],double eps,int n,double xx[]) {int i;double f1,f2,*x[2],ff,q,w;for(i=0;i<2;i++)x[i]=(double*)malloc(n*sizeof(double));for(i=0;i<n;i++){*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);*(x[1]+i)=a[i]+0.382*(b[i]-a[i]);}f1=objf(x[0]);f2=objf(x[1]);do{if(f1>f2){for(i=0;i<n;i++){b[i]=*(x[0]+i);*(x[0]+i)=*(x[1]+i);}f1=f2;for(i=0;i<n;i++)*(x[1]+i)=a[i]+0.382*(b[i]-a[i]);f2=objf(x[1]);}else{for(i=0;i<n;i++){a[i]=*(x[1]+i);*(x[1]+i)=*(x[0]+i);}f2=f1;for(i=0;i<n;i++)*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);f1=objf(x[0]);}q=0;for(i=0;i<n;i++)q=q+(b[i]-a[i])*(b[i]-a[i]);w=sqrt(q);}while(w>eps);for(i=0;i<n;i++)xx[i]=0.5*(a[i]+b[i]);ff=objf(xx);for(i=0;i<2;i++)free(x[i]);return(ff);}void jtf(double x0[],double h0,double s[],int n,double a[],double b[]) {int i;double*x[3],h,f1,f2,f3;for(i=0;i<3;i++)x[i]=(double*)malloc(n*sizeof(double));h=h0;for(i=0;i<n;i++)*(x[0]+i)=x0[i];f1=objf(x[0]);for(i=0;i<n;i++)*(x[1]+i)=*(x[0]+i)+h*s[i];f2=objf(x[1]);if(f2>=f1){h=-h0;for(i=0;i<n;i++)*(x[2]+i)=*(x[0]+i);f3=f1;for(i=0;i<n;i++){*(x[0]+i)=*(x[1]+i);*(x[1]+i)=*(x[2]+i);}f1=f2;f2=f3;}for(;;){h=2*h;for(i=0;i<n;i++)*(x[2]+i)=*(x[1]+i)+h*s[i];f3=objf(x[2]);if(f2<f3) break;else{for(i=0;i<n;i++){*(x[0]+i)=*(x[1]+i);*(x[1]+i)=*(x[2]+i);}f1=f2;f2=f3;}}if(h<0)for(i=0;i<n;i++){a[i]=*(x[2]+i);b[i]=*(x[0]+i);}elsefor(i=0;i<n;i++){a[i]=*(x[0]+i);b[i]=*(x[2]+i);}for(i=0;i<3;i++)free(x[i]);}double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[]) {double*a,*b,ff;a=(double*)malloc(n*sizeof(double));b=(double*)malloc(n*sizeof(double));jtf(x0,h0,s,n,a,b);ff=gold(a,b,epsg,n,x);free(a);free(b);return (ff);}double powell(double p[],double h0,double eps,double epsg,int n,double x[]) {int i,j,m;double*xx[4],*ss,*s;double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;ss=(double*)malloc(n*(n+1)*sizeof(double));s=(double*)malloc(n*sizeof(double));for(i=0;i<n;i++){for(j=0;j<=n;j++)*(ss+i*(n+1)+j)=0;*(ss+i*(n+1)+i)=1;}for(i=0;i<4;i++)xx[i]=(double*)malloc(n*sizeof(double));for(i=0;i<n;i++)*(xx[0]+i)=p[i];for(;;){for(i=0;i<n;i++){*(xx[1]+i)=*(xx[0]+i);x[i]=*(xx[1]+i);}f0=f1=objf(x);dlt=-1;for(j=0;j<n;j++){for(i=0;i<n;i++){*(xx[0]+i)=x[i];*(s+i)=*(ss+i*(n+1)+j);}f=oneoptim(xx[0],s,h0,epsg,n,x);df=f0-f;if(df>dlt){dlt=df;m=j;}}sdx=0;for(i=0;i<n;i++)sdx=sdx+fabs(x[i]-(*(xx[1]+i)));if(sdx<eps){free(ss);free(s);for(i=0;i<4;i++)free(xx[i]);return(f);}for(i=0;i<n;i++)*(xx[2]+i)=x[i];f2=f;for(i=0;i<n;i++){*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));x[i]=*(xx[3]+i);}fx=objf(x);f3=fx;q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);d=0.5*dlt*(f1-f3)*(f1-f3);if((f3<f1)||(q<d)){if(f2<=f3)for(i=0;i<n;i++)*(xx[0]+i)=*(xx[2]+i);elsefor(i=0;i<n;i++)*(xx[0]+i)=*(xx[3]+i);}else{for(i=0;i<n;i++){*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));*(s+i)=*(ss+(i+1)*(n+1));}f=oneoptim(xx[0],s,h0,epsg,n,x);for(i=0;i<n;i++)*(xx[0]+i)=x[i];for(j=m+1;j<=n;j++)for(i=0;i<n;i++)*(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);}}}void main(){double p[]={1,2};double ff,x[2];ff=powell(p,0.3,0.001,0.0001,2,x);printf("x[0]=%f,x[1]=%f,ff=%f\n",x[0],x[1],ff); getchar();}。
机械优化设计黄金分割法已知:F(x)=x4-4x3-6x2-16x+4,求极小值,极小值点,区间,迭代次数, 用进退法确定区间,用黄金分割法求极值。
#include <stdio.h>#include <math.h>#define e 0.001#define tt 0.01float f(double x){float y=pow(x,4)-4*pow(x,3)-6*pow(x,2)-16*x+4;return(y);}finding(float *p1,float*p2) {float x1=0,x2,x3,t,f1,f2,f3,h=tt; int n=0;x2=x1+h;f1=f(x1);f2=f(x2); if(f2>f1) {h=-h;t=x2;x2=x1;x1=t;} do { x3=x2+h;h=2*h;f3=f(x3);n=n+1;}while(f3<f2);if(x1>x3) {t=x1;x1=x3;x3=t;} *p1=x1;*p2=x3;return(n);}gold(float *p){float a,b,x1,x2,f1,f2; int n=0;finding(&a,&b);do{x1=a+0.382*(b-a);x2=a+0.618*(b-a);f1=f(x1);f2=f(x2);n=n+1;if(f1>f2) a=x1;else b=x2;}while((b-a)>e);*p=(x1+x2)/2;return(n); }main(){float a,b,x,min;int n1,n2; n1=finding(&a,&b);n2=gold(&x);min=f(x);printf("\n The area is %f to %f.",a,b); printf("\n The nunmber 1 is %d.",n1);printf("\n The min is %f and the result is %f.",x,min);printf("\n The nunmber 2 is %d.",n2)二插法已知:F(x1,x2)=4*x1-x2的平方-12;求极小值,极小值点,迭代次数,用复合形法求极值。
实验报告课程名称:机械优化设计实验项目:一维搜索(黄金分割)法上机实验专业班级: XXXXX级机械工程及自动化XX班学号: XXXXXXXXXX 姓名: XXXXXX 指导老师: XXXXXX 日期: 201X.12.12机械工程试验教学中心实验1 一维搜索(黄金分割)法实验报告实验日期 201X 年 12 月 11 日报告日期 201X 年 12 月 12 日班级 XXXXX级机自XXXX班姓名 XXXXXX 学号 XXXXXXXXXXXXXXX1、实验目的○1了解黄金分割法的基本原理;○2熟悉matlab程序使用方法;○3学习上机调试、运行所编写的程序。
2、黄金分割法原理该法适用于[a,b]区间上单谷函数极小值问题。
在搜索区间[a,b]内按照0.618比例加入两点α1,α2,并计算其函数值。
α1,α2将区间分成三段,然后利用区间消去法,通过比较函数值大小,删除其中一段,使搜索区间缩短,在保留区间进行同样处理,直到搜索区间缩小到指定精度为止。
3、编制MATLAB优化程序○1编写函数文件,并命名为fx.m保存,程序代码如下:function f=fx(w)%f=w^2-10*w+36;%f=w^4-5*w^3+4*w^2-6*w+60;%f=((w+1)^4)*((w-2)^2);注:上述“%”后面分别为要求解的三个方程,求解该方程式把相应方程式前面的“%”删除,点击保存,并运行下面的hjf.m文件,输入相应的初始步长h0、初始点x0、收敛法则epsilan的值○2编写进退法程序文件,命名为ab1.m保存,程序代码如下:function [a,b]=ab1(h0,x0)h=h0;x1=x0;f1=fx(x1);x2=x1+h;f2=fx(x2);if f2>f1h=-h;x3=x1;f3=f1;x1=x2;f1=f2;x2=x3;f2=f3;endx3=x2+h;f3=fx(x3);while f2>=f3x1=x2;f1=f2;x2=x3;f2=f3;x3=x2+2*h;f3=fx(x3);endif h<0a=x3;b=x1;elsea=x1;b=x3;end○3编写黄金分割法程序文件,命名为hjfgf.m 保存,程序代码如下: function hjfclearh1 = input('h0=?');x1=input('x0=?');epsilan=input('epsilan=?');[a,b]=ab1(h1,x1);x1=a+0.382*(b-a);f1=fx(x1);x2=a+0.618*(b-a);f2=fx(x2);while abs(b-a)>epsilanif f1>f2a=x1;x1=x2;f1=f2;x2=a+0.618*(b-a);f2=fx(x2);elseb=x2;x2=x1;f2=f1;x1=a+0.382*(b-a);f1=fx(x1);endendxm=(a+b)/2;['Optimal result:',blanks(3),'xm=[',...num2str(xm),']',blanks(6),'fm=',num2str(fx(xm))]4、实验结果()3610min )12+-=t t t f结果: *t = 5.0001 =*f 11 (h0=0.3、x0=0、epsilan=0.001)()60645min )2234+-+-=t t t t t f结果: *t = 3.2795 =*f 22.659 (h0=0.3、x0=0、epsilan=0.001)()()()2421min )3-+=t t t f 结果: *t =__-0.99989__ =*f __1.4253e-015__ (h0=0.3、x0=0、epsilan=0.001)。
#include"stdio.h"
#include"stdlib.h"
#include"math.h"
double objf(double x[])
{double ff;
ff=x[0]*x[0]+x[1]*x[1]-x[0]*x[1]-10*x[0]-4*x[1]+60;
return(ff);
}
double gold(double a[],double b[],double eps,int n,double xx[]) {int i;
double f1,f2,*x[2],ff,q,w;
for(i=0;i<2;i++)
x[i]=(double*)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
*(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
}
f1=objf(x[0]);
f2=objf(x[1]);
do
{if(f1>f2)
{for(i=0;i<n;i++)
{b[i]=*(x[0]+i);
*(x[0]+i)=*(x[1]+i);
}
f1=f2;
for(i=0;i<n;i++)
*(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
f2=objf(x[1]);
}
else
{for(i=0;i<n;i++)
{a[i]=*(x[1]+i);
*(x[1]+i)=*(x[0]+i);}
f2=f1;
for(i=0;i<n;i++)
*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
f1=objf(x[0]);
}
q=0;
for(i=0;i<n;i++)
q=q+(b[i]-a[i])*(b[i]-a[i]);
w=sqrt(q);
}while(w>eps);
for(i=0;i<n;i++)
xx[i]=0.5*(a[i]+b[i]);
ff=objf(xx);
for(i=0;i<2;i++)
free(x[i]);
return(ff);
}
void jtf(double x0[],double h0,double s[],int n,double a[],double b[]) {int i;
double*x[3],h,f1,f2,f3;
for(i=0;i<3;i++)
x[i]=(double*)malloc(n*sizeof(double));
h=h0;
for(i=0;i<n;i++)
*(x[0]+i)=x0[i];
f1=objf(x[0]);
for(i=0;i<n;i++)
*(x[1]+i)=*(x[0]+i)+h*s[i];
f2=objf(x[1]);
if(f2>=f1)
{h=-h0;
for(i=0;i<n;i++)
*(x[2]+i)=*(x[0]+i);
f3=f1;
for(i=0;i<n;i++)
{*(x[0]+i)=*(x[1]+i);
*(x[1]+i)=*(x[2]+i);
}
f1=f2;
f2=f3;
}
for(;;)
{h=2*h;
for(i=0;i<n;i++)
*(x[2]+i)=*(x[1]+i)+h*s[i];
f3=objf(x[2]);
if(f2<f3) break;
else
{for(i=0;i<n;i++)
{*(x[0]+i)=*(x[1]+i);
*(x[1]+i)=*(x[2]+i);
}
f1=f2;
f2=f3;
}
}
if(h<0)
for(i=0;i<n;i++)
{a[i]=*(x[2]+i);
b[i]=*(x[0]+i);
}
else
for(i=0;i<n;i++)
{a[i]=*(x[0]+i);
b[i]=*(x[2]+i);
}
for(i=0;i<3;i++)
free(x[i]);
}
double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[]) {double*a,*b,ff;
a=(double*)malloc(n*sizeof(double));
b=(double*)malloc(n*sizeof(double));
jtf(x0,h0,s,n,a,b);
ff=gold(a,b,epsg,n,x);
free(a);
free(b);
return (ff);
}
double powell(double p[],double h0,double eps,double epsg,int n,double x[]) {int i,j,m;
double*xx[4],*ss,*s;
double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;
ss=(double*)malloc(n*(n+1)*sizeof(double));
s=(double*)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{for(j=0;j<=n;j++)
*(ss+i*(n+1)+j)=0;
*(ss+i*(n+1)+i)=1;
}
for(i=0;i<4;i++)
xx[i]=(double*)malloc(n*sizeof(double));
for(i=0;i<n;i++)
*(xx[0]+i)=p[i];
for(;;)
{for(i=0;i<n;i++)
{*(xx[1]+i)=*(xx[0]+i);
x[i]=*(xx[1]+i);
}
f0=f1=objf(x);
dlt=-1;
for(j=0;j<n;j++)
{for(i=0;i<n;i++)
{*(xx[0]+i)=x[i];
*(s+i)=*(ss+i*(n+1)+j);
}
f=oneoptim(xx[0],s,h0,epsg,n,x);
df=f0-f;
if(df>dlt)
{dlt=df;
m=j;
}
}
sdx=0;
for(i=0;i<n;i++)
sdx=sdx+fabs(x[i]-(*(xx[1]+i)));
if(sdx<eps)
{free(ss);
free(s);
for(i=0;i<4;i++)
free(xx[i]);
return(f);
}
for(i=0;i<n;i++)
*(xx[2]+i)=x[i];
f2=f;
for(i=0;i<n;i++)
{*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));
x[i]=*(xx[3]+i);
}
fx=objf(x);
f3=fx;
q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);
d=0.5*dlt*(f1-f3)*(f1-f3);
if((f3<f1)||(q<d))
{if(f2<=f3)
for(i=0;i<n;i++)
*(xx[0]+i)=*(xx[2]+i);
else
for(i=0;i<n;i++)
*(xx[0]+i)=*(xx[3]+i);
}
else
{for(i=0;i<n;i++)
{*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));
*(s+i)=*(ss+(i+1)*(n+1));
}
f=oneoptim(xx[0],s,h0,epsg,n,x);
for(i=0;i<n;i++)
*(xx[0]+i)=x[i];
for(j=m+1;j<=n;j++)
for(i=0;i<n;i++)
*(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);
}
}
}
void main()
{double p[]={1,2};
double ff,x[2];
ff=powell(p,0.3,0.001,0.0001,2,x);
printf("x[0]=%f,x[1]=%f,ff=%f\n",x[0],x[1],ff); getchar();
}。