Kotaro Irisawa / Mbed 2 deprecated brent_only

2021.12.22

Dependencies:   mbed

main.cpp@1:5ec4a7674069, 2021-12-23 (annotated)

Committer:
Kotttaro
Date:
Thu Dec 23 06:57:42 2021 +0000
Revision:
1:5ec4a7674069
Parent:
0:904ead7b9c3a 
2021.12.23,15:57 succeed

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 1:5ec4a7674069 4 #define GOLD 1.618034
Kotttaro 1:5ec4a7674069 5 #define GLIMIT 100.0
Kotttaro 1:5ec4a7674069 6 #define TINY 1.0e-20
Kotttaro 0:904ead7b9c3a 7 #define SHIFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
Kotttaro 0:904ead7b9c3a 8 #define ZEPS 1.0e-10
Kotttaro 0:904ead7b9c3a 9 Serial pc(USBTX,USBRX);
Kotttaro 0:904ead7b9c3a 10 double f(double x);
Kotttaro 0:904ead7b9c3a 11 double brent(double min,double mid,double max,double tol);
Kotttaro 0:904ead7b9c3a 12 double SIGN(double x,double y);
Kotttaro 1:5ec4a7674069 13 void mnbrak();//引数は全てグローバルで指定するとする
Kotttaro 1:5ec4a7674069 14 double FMAX(double x,double y);
Kotttaro 0:904ead7b9c3a 15
Kotttaro 1:5ec4a7674069 16 double ax=0.0,bx=1.0,cx=0.0;
Kotttaro 1:5ec4a7674069 17 double fa,fb,fc;
Kotttaro 1:5ec4a7674069 18
Kotttaro 1:5ec4a7674069 19 void mnbrak()
Kotttaro 0:904ead7b9c3a 20 {
Kotttaro 1:5ec4a7674069 21 double ulim,u,r,q,fu,dum,fa,fb,fc;
Kotttaro 1:5ec4a7674069 22
Kotttaro 1:5ec4a7674069 23 fa=f(ax);
Kotttaro 1:5ec4a7674069 24 fb=f(bx);
Kotttaro 1:5ec4a7674069 25 pc.printf("fa=%lf,fb=%lf\r\n",fa,fb);
Kotttaro 1:5ec4a7674069 26 if(fb>fa)
Kotttaro 1:5ec4a7674069 27 {
Kotttaro 1:5ec4a7674069 28 SHIFT(dum,ax,bx,dum);
Kotttaro 1:5ec4a7674069 29 SHIFT(dum,fb,fa,dum);
Kotttaro 1:5ec4a7674069 30 }
Kotttaro 1:5ec4a7674069 31 cx=bx+GOLD*(bx-ax);
Kotttaro 1:5ec4a7674069 32 fc=f(cx);
Kotttaro 1:5ec4a7674069 33 while (fb>fc)
Kotttaro 1:5ec4a7674069 34 {
Kotttaro 1:5ec4a7674069 35 r=(bx-ax)*(fb-fc);
Kotttaro 1:5ec4a7674069 36 q=(bx-cx)*(fb-fa);
Kotttaro 1:5ec4a7674069 37 u=bx-((bx-cx)*q-(bx-cx)*r)/
Kotttaro 1:5ec4a7674069 38 (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));
Kotttaro 1:5ec4a7674069 39 ulim=bx+GLIMIT*(cx-bx);
Kotttaro 1:5ec4a7674069 40
Kotttaro 1:5ec4a7674069 41 if((bx-u)*(u-cx)>0.0)
Kotttaro 1:5ec4a7674069 42 {
Kotttaro 1:5ec4a7674069 43 fu=f(u);
Kotttaro 1:5ec4a7674069 44 if(fu<fc)
Kotttaro 1:5ec4a7674069 45 {
Kotttaro 1:5ec4a7674069 46 ax=bx;
Kotttaro 1:5ec4a7674069 47 bx=u;
Kotttaro 1:5ec4a7674069 48 fa=fb;
Kotttaro 1:5ec4a7674069 49 fb=fu;
Kotttaro 1:5ec4a7674069 50 return;
Kotttaro 1:5ec4a7674069 51 }
Kotttaro 1:5ec4a7674069 52 else if(fu>fb)
Kotttaro 1:5ec4a7674069 53 {
Kotttaro 1:5ec4a7674069 54 cx=u;
Kotttaro 1:5ec4a7674069 55 fc=fu;
Kotttaro 1:5ec4a7674069 56 return;
Kotttaro 1:5ec4a7674069 57 }
Kotttaro 1:5ec4a7674069 58 u=cx*+GOLD*(cx-bx);
Kotttaro 1:5ec4a7674069 59 fu=f(u);
Kotttaro 1:5ec4a7674069 60 }
Kotttaro 1:5ec4a7674069 61 else if((cx-u)*(u-ulim)>0.0)
Kotttaro 1:5ec4a7674069 62 {
Kotttaro 1:5ec4a7674069 63 fu=f(u);
Kotttaro 1:5ec4a7674069 64 if(fu<fc)
Kotttaro 1:5ec4a7674069 65 {
Kotttaro 1:5ec4a7674069 66 SHIFT(bx,cx,u,cx+GOLD*(cx-bx));
Kotttaro 1:5ec4a7674069 67 SHIFT(fb,fc,fu,f(u));
Kotttaro 1:5ec4a7674069 68 }
Kotttaro 1:5ec4a7674069 69
Kotttaro 1:5ec4a7674069 70 }
Kotttaro 1:5ec4a7674069 71 else if((u-ulim)*(ulim-cx)>=0.0)
Kotttaro 1:5ec4a7674069 72 {
Kotttaro 1:5ec4a7674069 73 u=ulim;
Kotttaro 1:5ec4a7674069 74 fu=f(u);
Kotttaro 1:5ec4a7674069 75 }
Kotttaro 1:5ec4a7674069 76 else
Kotttaro 1:5ec4a7674069 77 {
Kotttaro 1:5ec4a7674069 78 u=cx+GOLD*(cx-bx);
Kotttaro 1:5ec4a7674069 79 fu=f(u);
Kotttaro 1:5ec4a7674069 80 }
Kotttaro 1:5ec4a7674069 81 SHIFT(ax,bx,cx,u);
Kotttaro 1:5ec4a7674069 82 SHIFT(fa,fb,fc,fu);
Kotttaro 1:5ec4a7674069 83 }
Kotttaro 1:5ec4a7674069 84 }
Kotttaro 1:5ec4a7674069 85 double brent(double ax,double bx,double cx,double tol)
Kotttaro 1:5ec4a7674069 86 {
Kotttaro 1:5ec4a7674069 87 pc.printf("bernt start");
Kotttaro 0:904ead7b9c3a 88 int iter;
Kotttaro 0:904ead7b9c3a 89 double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm,xmin;
Kotttaro 0:904ead7b9c3a 90 double e=0.0;
Kotttaro 1:5ec4a7674069 91 a=(ax <cx ? ax : cx);
Kotttaro 1:5ec4a7674069 92 b=(ax >cx ? ax : cx);
Kotttaro 1:5ec4a7674069 93 x=w=v=bx;
Kotttaro 1:5ec4a7674069 94 fw=fv=fx=f(x);
Kotttaro 0:904ead7b9c3a 95 for(iter=1;iter<=ITMAX;iter++)
Kotttaro 0:904ead7b9c3a 96 {
Kotttaro 1:5ec4a7674069 97 //pc.printf("bernt loop");
Kotttaro 0:904ead7b9c3a 98 xm=0.5*(a+b);
Kotttaro 0:904ead7b9c3a 99 tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
Kotttaro 1:5ec4a7674069 100 pc.printf("x =%lf,w = %lf,u = %lf\n\r",x,w,u);
Kotttaro 0:904ead7b9c3a 101 if(fabs(x-xm)<=(tol2-0.5*(b-a)))
Kotttaro 0:904ead7b9c3a 102 {
Kotttaro 1:5ec4a7674069 103 //pc.printf("bernt out");
Kotttaro 0:904ead7b9c3a 104 xmin=x;
Kotttaro 1:5ec4a7674069 105 pc.printf("xmin=%lf\r\n",x);
Kotttaro 1:5ec4a7674069 106 pc.printf("fx=%lf\r\n",fx);
Kotttaro 0:904ead7b9c3a 107 return xmin;
Kotttaro 0:904ead7b9c3a 108 }
Kotttaro 0:904ead7b9c3a 109 if(fabs(e)>tol1)
Kotttaro 0:904ead7b9c3a 110 {
Kotttaro 0:904ead7b9c3a 111 r=(x-w)*(fx-fv);
Kotttaro 0:904ead7b9c3a 112 q=(x-v)*(fx-fw);
Kotttaro 0:904ead7b9c3a 113 p=(x-v)*q-(x-w)*r;
Kotttaro 0:904ead7b9c3a 114 q=2.0*(q-r);
Kotttaro 0:904ead7b9c3a 115 if(q>0.0)p=-p;
Kotttaro 0:904ead7b9c3a 116 q=fabs(q);
Kotttaro 0:904ead7b9c3a 117 etemp=e;
Kotttaro 0:904ead7b9c3a 118 e=d;
Kotttaro 0:904ead7b9c3a 119 if(fabs(p)>=fabs(0.5*q*etemp)||p<=q*(a-x)||p>=q*(b-x))
Kotttaro 0:904ead7b9c3a 120 { d=CGOLD*(e= (x>=xm ? a-x : b-x));}
Kotttaro 0:904ead7b9c3a 121 else
Kotttaro 0:904ead7b9c3a 122 {
Kotttaro 0:904ead7b9c3a 123 d=p/q;
Kotttaro 0:904ead7b9c3a 124 u=x+d;
Kotttaro 0:904ead7b9c3a 125 if(u-a < tol2 || b-u < tol2)
Kotttaro 0:904ead7b9c3a 126 {d=SIGN(tol1,xm-x);}
Kotttaro 0:904ead7b9c3a 127 }
Kotttaro 0:904ead7b9c3a 128 }
Kotttaro 0:904ead7b9c3a 129 else
Kotttaro 0:904ead7b9c3a 130 {
Kotttaro 0:904ead7b9c3a 131 d=CGOLD*(e= (x>=xm ? a-x : b-x));
Kotttaro 0:904ead7b9c3a 132 }
Kotttaro 0:904ead7b9c3a 133 u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
Kotttaro 1:5ec4a7674069 134 fu=f(u);///////注意
Kotttaro 0:904ead7b9c3a 135 if(fu <= fx)
Kotttaro 0:904ead7b9c3a 136 {
Kotttaro 0:904ead7b9c3a 137 if(u >= x)a=x; else b=x;
Kotttaro 0:904ead7b9c3a 138 SHIFT(v,w,x,u);
Kotttaro 0:904ead7b9c3a 139 SHIFT(fv,fw,fx,fu);
Kotttaro 0:904ead7b9c3a 140 }
Kotttaro 0:904ead7b9c3a 141 else{
Kotttaro 0:904ead7b9c3a 142 if(u < x){a=u;}
Kotttaro 0:904ead7b9c3a 143 else {b=u;}
Kotttaro 0:904ead7b9c3a 144 if(fu <= fw || w==x)
Kotttaro 0:904ead7b9c3a 145 {
Kotttaro 0:904ead7b9c3a 146 v=w;
Kotttaro 0:904ead7b9c3a 147 w=u;
Kotttaro 0:904ead7b9c3a 148 fv=fw;
Kotttaro 0:904ead7b9c3a 149 fw=fu;
Kotttaro 0:904ead7b9c3a 150 }
Kotttaro 0:904ead7b9c3a 151 else if (fu <= fv || v==x || v==w)
Kotttaro 0:904ead7b9c3a 152 {
Kotttaro 0:904ead7b9c3a 153 v=u;
Kotttaro 0:904ead7b9c3a 154 fv=fu;
Kotttaro 0:904ead7b9c3a 155 }
Kotttaro 0:904ead7b9c3a 156 }
Kotttaro 0:904ead7b9c3a 157
Kotttaro 0:904ead7b9c3a 158 }
Kotttaro 1:5ec4a7674069 159 pc.printf("xmin=%lf\r\n",x);
Kotttaro 1:5ec4a7674069 160 pc.printf("fx=%lf\r\n",fx);
Kotttaro 0:904ead7b9c3a 161 return xmin;
Kotttaro 0:904ead7b9c3a 162 }
Kotttaro 0:904ead7b9c3a 163 //極小値を求めたい関数を定義
Kotttaro 0:904ead7b9c3a 164 double f(double x){
Kotttaro 0:904ead7b9c3a 165 double x_return;
Kotttaro 1:5ec4a7674069 166 x_return=x*x*x-x;
Kotttaro 0:904ead7b9c3a 167 return x_return;
Kotttaro 0:904ead7b9c3a 168 }
Kotttaro 0:904ead7b9c3a 169 double SIGN(double x,double y)
Kotttaro 0:904ead7b9c3a 170 {
Kotttaro 0:904ead7b9c3a 171 double x_return;
Kotttaro 0:904ead7b9c3a 172 x_return=abs(x);
Kotttaro 0:904ead7b9c3a 173 if(y<0.0)x_return=-x_return;
Kotttaro 0:904ead7b9c3a 174 printf("f");
Kotttaro 0:904ead7b9c3a 175 return x_return;
Kotttaro 1:5ec4a7674069 176 }
Kotttaro 1:5ec4a7674069 177 double FMAX(double x, double y)
Kotttaro 1:5ec4a7674069 178 {
Kotttaro 1:5ec4a7674069 179 if(x>y){return x;}
Kotttaro 1:5ec4a7674069 180 if(y>x){return y;}
Kotttaro 1:5ec4a7674069 181 }
Kotttaro 0:904ead7b9c3a 182 int main() {
Kotttaro 1:5ec4a7674069 183 pc.baud(921600);
Kotttaro 1:5ec4a7674069 184 pc.printf("loop start\r\n");
Kotttaro 1:5ec4a7674069 185 mnbrak();
Kotttaro 1:5ec4a7674069 186 pc.printf("(a,b,c)=(%lf,%lf,%lf)\r\n",ax,bx,cx);
Kotttaro 1:5ec4a7674069 187 pc.printf("(fa,fb,fc)=(%lf,%lf,%lf)\r\n",fa,fb,fc);
Kotttaro 1:5ec4a7674069 188 brent(ax,bx,cx,0.000001);
Kotttaro 0:904ead7b9c3a 189 return 0;
Kotttaro 0:904ead7b9c3a 190 }