solve .

“msharpmath”, The Simple is the Best

[100] 005 Tutorial solve           root and minimum finding

Click the above button for full-text tutorial of   ‘solve’ ,  ‘linsys’ and ‘minfind’.

//==================================================================================
//   Umbrella  ‘solve’    for coupled nonlinear equations.
//   Umbrella  ‘linsys’   for pseudo-linear system of nonlinear equations
//   Umbrella  ‘minfind’   for a local minimum
//==================================================================================

//==================================================================================
//   syntax
//———————————————————————————-
//  solve .x { a } ( <<opt>>, f(x) = g(x) )  // default method
//  solve .x         ( <<opt>>, f(x) = g(x) )  // default method with a = 0.41

The initial guess x = a by {a} can be omitted and a default value a = 0.41 is given instead.

//  solve .x ( <<opt>>, f(x) = g(x) ) .bisect( a,b )              // bisection method
//  solve .x ( <<opt>>, f(x) = g(x) ) .secant( a,b )              // secant method
//  solve .x ( <<opt>>, f(x) = g(x) ) .muller( a,b,c )            // Muller method
//  solve .x ( <<opt>>, f(x) = g(x) ) .march ( a,b,dx )           // marching method
//  solve .x ( <<opt>>, f(x) = g(x) ) .span[n=100,g=1](a,b) // bisection over span

//  solve .x { a } ( <<opt>>, f(x) = g(x) ) .newton  ( f’(x)-g’(x) )  // Newton method
//  solve .x { a } ( <<opt>>, f(x) = g(x) ) .multiple( f’(x)-g’(x) )  // multiple-roots
//———————————————————————————-
// Spokes

//  .peep                show all the iteration procedures
//  .peeptail            show the final stage of the iteration
//  .return( … )        return manipulated data
//  .maxiter(1000)       maximum number of iteration
//  .tol/abstol(1.e-7)   absolute tolerance. Either tol or abstol is equivalent
//  .reltol(1.e-7)       relative tolerance

#> solve .x { 0.5 } ( exp(x)-3 = sin(x) );
ans =       1.3818344

#> solve .x { 0 }   ( exp(x)-3 = sin(x) );
The secant method fails due to near-zero slope
ans =          1.#INF

%> Spoke return============================================

#> solve .x { 0.5 } ( exp(x)-3 = sin(x) ) .return( x/100 );
ans =     0.013818344

#> solve .x { 0.5 } ( exp(x)-3 = sin(x) ) .return( x+1, x, sin(x) );
ans =
[       2.38183 ]
[       1.38183 ]
[        0.9822 ]

#> solve .x ( exp(x)-3 = sin(x) ) .bisect( 1,2 ).peeptail;
——— bisection method ——–
iter              x         f(x)
23       1.3818344   7.724e-008
ans =       1.3818344

#> solve .x ( exp(x)-3-sin(x) = 0 ) .secant( 1,2 ).peep;
——— secant method ———–
iter              x         f(x)
0               1       -1.123
0               2         3.48
1       1.2440152      -0.4776
2       1.4245115       0.1665
3        1.377849     -0.01508
4       1.3817247   -0.0004161
5       1.3818347   1.087e-006
6       1.3818344  -7.796e-011
ans =       1.3818344

#> solve .x ( exp(x)-3 = sin(x) ) .march( 1,2, 0.1 ).peep;
——— marching method ———–
step              x         f(x)
0               1       -1.123
1             1.1       -0.887
2             1.2      -0.6119
3             1.3      -0.2943
4             1.4      0.06975
ans =       1.3818343

//==================================================================================
//  Collection of Roots
//==================================================================================
#> solve .x ( sin(x) = 0.5 ) .span(0,5*pi) /pi*6;
ans =
[             1  5.00549e-008 ]
[             5 -1.06018e-016 ]
[            13 -6.00659e-007 ]
[            17 -6.00659e-007 ]
[            25  1.06018e-015 ]
[            29 -1.20132e-006 ]

#> solve .x ( .J_0(3*x)-2*sin(4*x)^2 = 0 ) .span[101](0,5);
ans =
[      0.186725  7.04574e-008 ]
[       0.72046 -4.88473e-008 ]
[      0.798234 -8.11431e-008 ]
[       2.25829  -8.5721e-007 ]
[       2.45261 -9.71382e-007 ]
[       4.63561 -1.21728e-006 ]
[       4.77437  1.29266e-006 ]

#> plot.x[1001](0,5) ( .J_0(3*x)-2*sin(4*x)^2, 0 );

%> From the plot, it is possible to have roots near x = 4. Therefore, more dense span is applied as follows (from 100 span to 1000 span)

#> solve .x ( .J_0(3*x)-2*sin(4*x)^2 = 0 ) .span[1001](0,5);
ans =
[      0.186725 -5.51991e-008 ]
[       0.72046  3.84365e-008 ]
[      0.798234   9.7982e-009 ]
[       2.25829  2.61073e-007 ]
[       2.45261 -2.99433e-007 ]
[       3.93141 -1.25321e-008 ]
[       3.94436  3.00444e-008 ]
[       4.63561 -8.92462e-007 ]
[       4.77437 -3.62101e-007 ]

//==================================================================================
//  Complex Root
//==================================================================================

#> solve .z { 1+10! } ( exp(z)+1 = 0 ).peep / pi;
———— secant method for complex root ————
iter             real             imag          |f(z)|
0  (              1               10 )         1.956
0  (          1.011            10.03 )         2.008
1  (     0.65476666        9.9056557 )         1.136
2  (     0.42414361        9.8163938 )        0.7145
3  (     0.24343609        9.7135918 )        0.4262
4  (     0.12477307        9.6178385 )        0.2444
5  (    0.057967283        9.5412233 )        0.1338
6  (    0.025610122        9.4894157 )       0.07041
7  (    0.011366113        9.4587996 )       0.03607
8  (   0.0052156209        9.4421955 )       0.01823
9  (   0.0024756637        9.4335815 )      0.009156
10  (   0.0012028686        9.4292021 )      0.004587
11  (  0.00059245921        9.4269954 )      0.002296
12  (  0.00029395732         9.425888 )      0.001148
13  (  0.00014640709        9.4253333 )     0.0005744
14  ( 7.3060216e-005        9.4250557 )     0.0002872
15  ( 3.6494222e-005        9.4249169 )     0.0001436
16  ( 1.8238133e-005        9.4248474 )    7.181e-005
17  ( 9.1168212e-006        9.4248127 )    3.591e-005
18  ( 4.5578491e-006        9.4247953 )    1.795e-005
19  ( 2.2787842e-006        9.4247866 )    8.977e-006
20  (  1.139357e-006        9.4247823 )    4.488e-006
21  ( 5.6966973e-007        9.4247801 )    2.244e-006
22  ( 2.8483267e-007         9.424779 )    1.122e-006
23  ( 1.4241579e-007        9.4247785 )     5.61e-007
24  ( 7.1207756e-008        9.4247782 )    2.805e-007
25  ( 3.5603844e-008        9.4247781 )    1.403e-007
26  ( 1.7801913e-008         9.424778 )    7.013e-008
ans =  5.66653e-009 + 3!

#> solve .z { 1+20! } ( exp(z)+1 = 0 ) / pi;
ans =  -1.73861e-008 + 7!

//==================================================================================
//  Coupled Nonlinear Equations
//==================================================================================
// syntax
//———————————————————————————-
//  solve .x.y.z … { xo,yo,zo, … } ( <<opt>>, f = 0, g = 0, h = 0, … )
//  solve [ x=xo, y=yo, z=zo ] ( <<opt>>, f = 0, g = 0, h = 0, … )
//———————————————————————————-
// Spokes
//———————————————————————————-
//  .relax(0.7)          relaxation factor to improve convergence
//  .broyden             use Broyden method instead of the default method

#> solve .x.y { 1,1 } ( x*x*y+y^3 = 39, x*x+3*x*y-y*y = 13 ) .peep;
— Newton method with relaxation (default 0.7)
iter      residue           |dx|
0         38.33
1         380.5          9.168
2           169          2.416
3         72.73          1.706
4         29.68          1.135
5         11.16         0.6685
6         3.811         0.3258
7          1.21         0.1285
8        0.3702        0.04331
9        0.1117        0.01353
10       0.03359       0.004113
11       0.01008       0.001239
12      0.003025      0.0003721
13     0.0009077      0.0001117
14     0.0002723      3.35e-005
15     8.17e-005     1.005e-005
16    2.451e-005     3.015e-006
17    7.353e-006     9.047e-007
18    2.206e-006     2.714e-007
19    6.618e-007     8.142e-008
ans =
[             2 ]
[             3 ]

#> solve [ x=1,y=1 ] ( x*x*y+y^3 = 39, x*x+3*x*y-y*y = 13 ) ;  x; y;
ans =
[             2 ]
[             3 ]
x =               2
y =               3

#> solve .F.u { 0.1,1 } (
u*u*(1+750*F^2) = 300,
1 = F*( 2*log10(4.e4*fabs(u*F))-0.8 )
).peep
.return( F*F, u, F, u*u*(1+750*F^2)-300, -1+F*(2*log10(4.e4*u*F)-0.8) ) ;
— Newton method with relaxation (default 0.7)
iter      residue           |dx|
0         291.5
1          32.9          18.68
2          5531          25.29
3          1401          25.55
4         468.8          7.749
5         159.6            3.2
6         51.79          1.301
7         16.06         0.4674
8         4.874         0.1512
9         1.468        0.04654
10        0.4407        0.01407
11        0.1323       0.004233
12       0.03968       0.001271
13       0.01191      0.0003813
14      0.003572      0.0001144
15      0.001071     3.432e-005
16     0.0003214      1.03e-005
17    9.643e-005     3.089e-006
18    2.893e-005     9.267e-007
19    8.679e-006      2.78e-007
20    2.604e-006      8.34e-008
ans =
[     0.0157394 ]
[       4.84036 ]
[      0.125457 ]
[  2.60368e-006 ]
[ -4.34674e-010 ]

#> solve .x.y ( 3*x+ 6*y = 25+5*y, 2*x-4*y = 2-3*x ) .relax(1).peep;
— Newton method with relaxation (default 0.7)
iter      residue           |dx|
0         23.41
1    1.473e-009          8.634
2    7.944e-015     5.728e-010
ans =
[             6 ]
[             7 ]

//==================================================================================
//  Pseudo-Linear System
//==================================================================================
// syntax
//  linsys .x.y. …  ( <<opt>>, f(x,y,…),g(x,y,…), … )
//  linsys .x.y. … { a,b, … } ( <<opt>>, f(x,y,…),g(x,y,…), … )
// linsys [ x=a, y=b, … ] ( <<opt>>, f(x,y,…),g(x,y,…), … )
// x=a; y=b; …  linsys [ x, y, … ] ( <<opt>>, f(x,y,…),g(x,y,…), … )

#> x = y = 1;;

#> for.i(1,4) linsys [x,y] ( x*x+y*y = 4, 2*x-y = 1 );
ans =  [       1.33333       1.66666 ]
ans =  [       1.27381       1.54762 ]
ans =  [       1.27178       1.54356 ]
ans =  [       1.27178       1.54356 ]

//==================================================================================
//  Mininum Search
//==================================================================================
//———————————————————————————-
// syntax
//———————————————————————————-
//  minfind .x1.x2 … .xn { c1,c2, … , cn } (  <<opt>>, F(x1,x2, … ,xn) )
//  minfind [ x1=c1, x2=c2, … , xn=cn ]     (  <<opt>>, F(x1,x2, … ,xn) )
//  minfind .x1.x2 … .xn { c1,c2, … , cn } (  <<opt>>, f1,f2, … , fn )
//———————————————————————————-
//  Spokes
//———————————————————————————-
//  .tol/abstol/reltol(1.e-5)    tolerance to terminate iteration
//  .maxiter(200)                maximum iteration permitted
//  .maxstep(1)                  maximum step size toward negative gradient direction
//  .peep                        show the iteration procedure
//  .peeptail                    show the final stage of iteration
//  .golden                      use the golden-section method (this is default)
//  .quadratic                   use the quadratic method
//  .fletcher                    use the Fletcher-Powell method

#> minfind .x.y {1,1}  ( 5*x*x -6*x*y +2*y*y -2*y +4 );  // x,y not defined
ans =
[       2.99034 ]
[       4.98455 ]
[     -0.999933 ]

#> minfind [ x=1,y=1 ] ( 5*x*x -6*x*y +2*y*y -2*y +4 ); x; y;
ans =
[       2.99034 ]
[       4.98455 ]
[     -0.999933 ]
x =       2.9903372
y =       4.9845459

#> minfind.x.y {1,1} (  x*x+y*y = 4, 2*x-y = 1 )
.return( x,y, x*x+y*y-4, 2*x-y-1 );
ans =
[       1.27196 ]
[       1.54336 ]
[  6.38006e-006 ]
[   0.000560276 ]

Comments are closed.