Kotaro Irisawa / Mbed 2 deprecated brent_only

2021.12.22

Dependencies:   mbed

main.cpp@0:904ead7b9c3a, 2021-12-22 (annotated)

Committer:
Kotttaro
Date:
Wed Dec 22 12:42:04 2021 +0000
Revision:
0:904ead7b9c3a
Child:
1:5ec4a7674069
2021.12.22;

Who changed what in which revision?

UserRevisionLine numberNew contents of line
Kotttaro 0:904ead7b9c3a 1 #include "mbed.h"
Kotttaro 0:904ead7b9c3a 2 #define ITMAX 100
Kotttaro 0:904ead7b9c3a 3 #define CGOLD 0.3819660
Kotttaro 0:904ead7b9c3a 4 #define SHIFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
Kotttaro 0:904ead7b9c3a 5 #define ZEPS 1.0e-10
Kotttaro 0:904ead7b9c3a 6 Serial pc(USBTX,USBRX);
Kotttaro 0:904ead7b9c3a 7 double f(double x);
Kotttaro 0:904ead7b9c3a 8 double brent(double min,double mid,double max,double tol);
Kotttaro 0:904ead7b9c3a 9 double SIGN(double x,double y);
Kotttaro 0:904ead7b9c3a 10
Kotttaro 0:904ead7b9c3a 11 double brent(double min,double mid,double max,double tol)
Kotttaro 0:904ead7b9c3a 12 {
Kotttaro 0:904ead7b9c3a 13 int iter;
Kotttaro 0:904ead7b9c3a 14 double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm,xmin;
Kotttaro 0:904ead7b9c3a 15 double e=0.0;
Kotttaro 0:904ead7b9c3a 16 a=(min < max ? min : max);
Kotttaro 0:904ead7b9c3a 17 b=(min > max ? min : max);
Kotttaro 0:904ead7b9c3a 18 x=w=v=mid;
Kotttaro 0:904ead7b9c3a 19 fw=fv=fu=f(x);
Kotttaro 0:904ead7b9c3a 20 for(iter=1;iter<=ITMAX;iter++)
Kotttaro 0:904ead7b9c3a 21 {
Kotttaro 0:904ead7b9c3a 22 xm=0.5*(a+b);
Kotttaro 0:904ead7b9c3a 23 tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
Kotttaro 0:904ead7b9c3a 24 if(fabs(x-xm)<=(tol2-0.5*(b-a)))
Kotttaro 0:904ead7b9c3a 25 {
Kotttaro 0:904ead7b9c3a 26 xmin=x;
Kotttaro 0:904ead7b9c3a 27 return xmin;
Kotttaro 0:904ead7b9c3a 28 }
Kotttaro 0:904ead7b9c3a 29 if(fabs(e)>tol1)
Kotttaro 0:904ead7b9c3a 30 {
Kotttaro 0:904ead7b9c3a 31 r=(x-w)*(fx-fv);
Kotttaro 0:904ead7b9c3a 32 q=(x-v)*(fx-fw);
Kotttaro 0:904ead7b9c3a 33 p=(x-v)*q-(x-w)*r;
Kotttaro 0:904ead7b9c3a 34 q=2.0*(q-r);
Kotttaro 0:904ead7b9c3a 35 if(q>0.0)p=-p;
Kotttaro 0:904ead7b9c3a 36 q=fabs(q);
Kotttaro 0:904ead7b9c3a 37 etemp=e;
Kotttaro 0:904ead7b9c3a 38 e=d;
Kotttaro 0:904ead7b9c3a 39 if(fabs(p)>=fabs(0.5*q*etemp)||p<=q*(a-x)||p>=q*(b-x))
Kotttaro 0:904ead7b9c3a 40 { d=CGOLD*(e= (x>=xm ? a-x : b-x));}
Kotttaro 0:904ead7b9c3a 41 else
Kotttaro 0:904ead7b9c3a 42 {
Kotttaro 0:904ead7b9c3a 43 d=p/q;
Kotttaro 0:904ead7b9c3a 44 u=x+d;
Kotttaro 0:904ead7b9c3a 45 if(u-a < tol2 || b-u < tol2)
Kotttaro 0:904ead7b9c3a 46 {d=SIGN(tol1,xm-x);}
Kotttaro 0:904ead7b9c3a 47 }
Kotttaro 0:904ead7b9c3a 48 }
Kotttaro 0:904ead7b9c3a 49 else
Kotttaro 0:904ead7b9c3a 50 {
Kotttaro 0:904ead7b9c3a 51 d=CGOLD*(e= (x>=xm ? a-x : b-x));
Kotttaro 0:904ead7b9c3a 52 }
Kotttaro 0:904ead7b9c3a 53 u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
Kotttaro 0:904ead7b9c3a 54 fu=f(x);
Kotttaro 0:904ead7b9c3a 55 if(fu <= fx)
Kotttaro 0:904ead7b9c3a 56 {
Kotttaro 0:904ead7b9c3a 57 if(u >= x)a=x; else b=x;
Kotttaro 0:904ead7b9c3a 58 SHIFT(v,w,x,u);
Kotttaro 0:904ead7b9c3a 59 SHIFT(fv,fw,fx,fu);
Kotttaro 0:904ead7b9c3a 60 }
Kotttaro 0:904ead7b9c3a 61 else{
Kotttaro 0:904ead7b9c3a 62 if(u < x){a=u;}
Kotttaro 0:904ead7b9c3a 63 else {b=u;}
Kotttaro 0:904ead7b9c3a 64 if(fu <= fw || w==x)
Kotttaro 0:904ead7b9c3a 65 {
Kotttaro 0:904ead7b9c3a 66 v=w;
Kotttaro 0:904ead7b9c3a 67 w=u;
Kotttaro 0:904ead7b9c3a 68 fv=fw;
Kotttaro 0:904ead7b9c3a 69 fw=fu;
Kotttaro 0:904ead7b9c3a 70 }
Kotttaro 0:904ead7b9c3a 71 else if (fu <= fv || v==x || v==w)
Kotttaro 0:904ead7b9c3a 72 {
Kotttaro 0:904ead7b9c3a 73 v=u;
Kotttaro 0:904ead7b9c3a 74 fv=fu;
Kotttaro 0:904ead7b9c3a 75 }
Kotttaro 0:904ead7b9c3a 76 }
Kotttaro 0:904ead7b9c3a 77
Kotttaro 0:904ead7b9c3a 78 }
Kotttaro 0:904ead7b9c3a 79 return xmin;
Kotttaro 0:904ead7b9c3a 80 }
Kotttaro 0:904ead7b9c3a 81 //極小値を求めたい関数を定義
Kotttaro 0:904ead7b9c3a 82 double f(double x){
Kotttaro 0:904ead7b9c3a 83 double x_return;
Kotttaro 0:904ead7b9c3a 84 x_return=x*x;
Kotttaro 0:904ead7b9c3a 85 return x_return;
Kotttaro 0:904ead7b9c3a 86 }
Kotttaro 0:904ead7b9c3a 87 double SIGN(double x,double y)
Kotttaro 0:904ead7b9c3a 88 {
Kotttaro 0:904ead7b9c3a 89 double x_return;
Kotttaro 0:904ead7b9c3a 90 x_return=abs(x);
Kotttaro 0:904ead7b9c3a 91 if(y<0.0)x_return=-x_return;
Kotttaro 0:904ead7b9c3a 92 printf("f");
Kotttaro 0:904ead7b9c3a 93 return x_return;
Kotttaro 0:904ead7b9c3a 94 }
Kotttaro 0:904ead7b9c3a 95 int main() {
Kotttaro 0:904ead7b9c3a 96 double x;
Kotttaro 0:904ead7b9c3a 97 x= brent(-100.0,10.0,100.0,0.0000000001);
Kotttaro 0:904ead7b9c3a 98 printf("%lf\r\n",x);
Kotttaro 0:904ead7b9c3a 99 return 0;
Kotttaro 0:904ead7b9c3a 100 }