Kotaro Irisawa
/
brent_only
2021.12.22
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?
User | Revision | Line number | New 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 | } |