Revision 1:5ec4a7674069, committed 2021-12-23

Comitter:

Kotttaro

Date:

Thu Dec 23 06:57:42 2021 +0000

Parent:

0:904ead7b9c3a

Commit message:

2021.12.23,15:57 succeed

Changed in this revision

--- a/main.cpp	Wed Dec 22 12:42:04 2021 +0000
+++ b/main.cpp	Thu Dec 23 06:57:42 2021 +0000
@@ -1,29 +1,109 @@
 #include "mbed.h"
 #define ITMAX 100
 #define CGOLD 0.3819660
+#define GOLD 1.618034
+#define GLIMIT 100.0
+#define TINY 1.0e-20
 #define SHIFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
 #define ZEPS 1.0e-10
 Serial pc(USBTX,USBRX);
 double f(double x);
 double brent(double min,double mid,double max,double tol);
 double SIGN(double x,double y);
+void mnbrak();//引数は全てグローバルで指定するとする
+double FMAX(double x,double y);
 
-double brent(double min,double mid,double max,double tol)
+double ax=0.0,bx=1.0,cx=0.0;
+double fa,fb,fc;
+
+void mnbrak()
 {
+    double ulim,u,r,q,fu,dum,fa,fb,fc;
+    
+    fa=f(ax);
+    fb=f(bx);
+    pc.printf("fa=%lf,fb=%lf\r\n",fa,fb);
+    if(fb>fa)
+    {
+        SHIFT(dum,ax,bx,dum);
+        SHIFT(dum,fb,fa,dum);
+        }
+    cx=bx+GOLD*(bx-ax);
+    fc=f(cx);
+    while (fb>fc)
+    {
+        r=(bx-ax)*(fb-fc);
+        q=(bx-cx)*(fb-fa);
+        u=bx-((bx-cx)*q-(bx-cx)*r)/
+        (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));
+        ulim=bx+GLIMIT*(cx-bx);
+        
+        if((bx-u)*(u-cx)>0.0)
+        {
+            fu=f(u);
+            if(fu<fc)
+            {
+                ax=bx;
+                bx=u;
+                fa=fb;
+                fb=fu;
+                return;
+                }
+            else if(fu>fb)
+            {
+                cx=u;
+                fc=fu;
+                return;
+                }
+             u=cx*+GOLD*(cx-bx);
+             fu=f(u);       
+            }
+        else if((cx-u)*(u-ulim)>0.0)
+        {
+            fu=f(u);
+            if(fu<fc)
+            {
+                SHIFT(bx,cx,u,cx+GOLD*(cx-bx));
+                SHIFT(fb,fc,fu,f(u));
+                }
+            
+            }
+        else if((u-ulim)*(ulim-cx)>=0.0)
+        {
+            u=ulim;
+            fu=f(u);
+            }
+        else
+        {
+            u=cx+GOLD*(cx-bx);
+            fu=f(u);
+            }
+        SHIFT(ax,bx,cx,u);
+        SHIFT(fa,fb,fc,fu);
+        }
+    }
+double brent(double ax,double bx,double cx,double tol)
+{
+    pc.printf("bernt start");
     int iter;
     double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm,xmin;
     double e=0.0;
-    a=(min < max ? min : max);
-    b=(min > max ? min : max);    
-    x=w=v=mid;
-    fw=fv=fu=f(x);
+    a=(ax <cx ? ax : cx);
+    b=(ax >cx ? ax : cx);
+    x=w=v=bx;
+    fw=fv=fx=f(x);
     for(iter=1;iter<=ITMAX;iter++)
     {
+        //pc.printf("bernt loop");
         xm=0.5*(a+b);
         tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
+        pc.printf("x =%lf,w = %lf,u = %lf\n\r",x,w,u);
         if(fabs(x-xm)<=(tol2-0.5*(b-a)))
         {
+            //pc.printf("bernt out");
             xmin=x;
+             pc.printf("xmin=%lf\r\n",x);
+             pc.printf("fx=%lf\r\n",fx);
             return xmin;
         }
         if(fabs(e)>tol1)
@@ -51,7 +131,7 @@
             d=CGOLD*(e= (x>=xm ? a-x : b-x));
             }    
         u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
-        fu=f(x);
+        fu=f(u);///////注意
          if(fu <= fx)
          {
             if(u >= x)a=x; else b=x;
@@ -76,12 +156,14 @@
         }
         
         }
+        pc.printf("xmin=%lf\r\n",x);
+        pc.printf("fx=%lf\r\n",fx);
     return xmin;
     }
 //極小値を求めたい関数を定義    
 double f(double x){
     double x_return;
-    x_return=x*x;
+    x_return=x*x*x-x;
     return x_return;
     }   
 double SIGN(double x,double y)
@@ -91,10 +173,18 @@
     if(y<0.0)x_return=-x_return;
     printf("f");
     return x_return; 
-    }     
+    } 
+double FMAX(double x, double y)
+{
+    if(x>y){return x;}
+    if(y>x){return y;}
+    }        
 int main() {
-   double x;
-   x= brent(-100.0,10.0,100.0,0.0000000001);
-   printf("%lf\r\n",x);
+   pc.baud(921600);
+   pc.printf("loop start\r\n");
+   mnbrak();
+   pc.printf("(a,b,c)=(%lf,%lf,%lf)\r\n",ax,bx,cx);
+   pc.printf("(fa,fb,fc)=(%lf,%lf,%lf)\r\n",fa,fb,fc);
+   brent(ax,bx,cx,0.000001);
    return 0;
 }