Small Solver Program

02202019, 05:31 AM
Post: #21




RE: Small Solver Program
(02192019 08:39 AM)Csaba Tizedes Wrote: I have another idea: what if, if we just simple walking along on the function's curve until the sign is changing, then we turning back, decreasing the stepsize and walking again. If the sign of the function changed, we are turning back again, decreasing the size of the steps and so on... If the stepsize is less than a small value, we found the root. Or if we start with two points whose function values have opposite signs and then evaluate the function at the midpoint. Based on its sign the next step is just a jump to the left or a jump to the right. Bang: we've just invented binary search. Cheers Thomas 

02202019, 07:22 AM
(This post was last modified: 02202019 08:39 AM by Thomas Klemm.)
Post: #22




RE: Small Solver Program
(02182019 10:22 PM)Thomas Klemm Wrote: Which I still do not understand in detail. (02192019 01:10 AM)Albert Chan Wrote: I think the code add slope data, like this: If we know two points \((x_1, y_1)\) and \((x_2, y_2)\) on a parabola \(y=ax^2+bx+c\) then we can calculate the slope \(y'=2ax+b\) at their midpoint \(\frac{x_1 + x_2}{2}\): \(\begin{matrix} y_1 = ax_1^2+bx_1+c \\ y_2 = ax_2^2+bx_2+c \end{matrix}\) \(\begin{align*} y_2y_1 &= a(x_2^2x_1^2)+b(x_2x_1) \\ &= a(x_2+x_1)(x_2x_1)+b(x_2x_1) \end{align*}\) \(\begin{align*} \frac{y_2y_1}{x_2x_1} &= a(x_2+x_1)+b \\ &= 2a\frac{x_2+x_1}{2}+b \end{align*}\) This is a bit surprising since a parabola is determined by three points, but the result is independent of the third point. With linear regression we try to minimise the distance of the vector \(y=\begin{bmatrix} y_1\\ \cdot\\ \cdot\\ \cdot\\ y_n \end{bmatrix}\) to the 3dimensional subspace spanned by \(x^2=\begin{bmatrix} x_1^2\\ \cdot\\ \cdot\\ \cdot\\ x_n^2 \end{bmatrix}\, , \, x=\begin{bmatrix} x_1\\ \cdot\\ \cdot\\ \cdot\\ x_n \end{bmatrix}\, , \, e=\begin{bmatrix} 1\\ \cdot\\ \cdot\\ \cdot\\ 1 \end{bmatrix}\) The solution is the orthogonal projection. Now we take Csaba's transformation: \(\Delta= \begin{bmatrix} \frac{1}{x_2x_1} & \frac{1}{x_2x_1} & 0 & 0 & \cdots & 0 \\ 0 & \frac{1}{x_3x_2} & \frac{1}{x_3x_2} & 0 & \cdots & 0 \\ \cdot & \cdot & \cdot & \cdot & & \cdot \\ \cdot & \cdot & \cdot & \cdot & & \cdot \\ \cdot & \cdot & \cdot & \cdot & & \cdot \\ 0 & 0 & 0 & \cdots & \frac{1}{x_nx_{n1}} & \frac{1}{x_nx_{n1}} \end{bmatrix}\) and apply it to all the vectors \(\Delta y=\begin{bmatrix} \frac{y_2y_1}{x_2x1} \\ \cdot\\ \cdot\\ \cdot\\ \frac{y_ny_{n1}}{x_nx_{n1}} \end{bmatrix} \) \(\Delta x^2=\begin{bmatrix} \frac{x_2^2x_1^2}{x_2x1} \\ \cdot\\ \cdot\\ \cdot\\ \frac{x_n^2x_{n1}^2}{x_nx_{n1}} \end{bmatrix}=\begin{bmatrix} x_2+x_1\\ \cdot\\ \cdot\\ \cdot\\ x_n+x_{n1} \end{bmatrix}\) \(\Delta x=\begin{bmatrix} \frac{x_2x_1}{x_2x1} \\ \cdot\\ \cdot\\ \cdot\\ \frac{x_nx_{n1}}{x_nx_{n1}} \end{bmatrix}=\begin{bmatrix} 1\\ \cdot\\ \cdot\\ \cdot\\ 1 \end{bmatrix}\) \(\Delta e=\begin{bmatrix} 0\\ \cdot\\ \cdot\\ \cdot\\ 0 \end{bmatrix}\) Since \(\Delta e = 0\) we realise that we project along the vector \(e\) into a n1 dimensional subspace. Here we try to minimise the distance of the vector \(\Delta y\) to the 2dimensional subspace spanned by \(\Delta x^2\) and \(\Delta x\). Again the solution is the orthogonal projection. But clearly the transformation \(\Delta\) isn't orthogonal. Thus these two solutions aren't mapped exactly. What reduces this error? Csaba suggested sorting the points by their xcoordinates. In case of his example this leads to an "isotropic" transformation \(\Delta\) since all entries of the matrix are the same except those that are 0 of course. Is that why it works so well? Is it enough to keep them in order, or do we also need the same distances? I could imagine that this leads to a completely wrong result if the transformation \(\Delta\) is a distortion. Cheers Thomas 

02242019, 09:21 AM
Post: #23




RE: Small Solver Program
(02182019 10:22 PM)Thomas Klemm Wrote: Quadratic fit without linear algebra (32SII)  NOT only for statistics lovers  LONG (02192019 08:39 AM)Csaba Tizedes Wrote: But as a programming question this was an interesting "game" for me. Since the regression is linear you can remove division by 2 at two places in your program. First when calculating the midpoints: Code: RCL X And then later on when calculating A: Code: LBL C #Calculating A, B and C We can rewrite \(y=m\frac{x_k+x_{k+1}}{2}+b\) as \(y=\frac{m}{2}(x_k+x_{k+1})+b\). But since \(A=\frac{m}{2}\) this becomes \(y=A(x_k+x_{k+1})+b\). Cheers Thomas 

02252019, 08:39 PM
Post: #24




RE: Small Solver Program
(02202019 05:31 AM)Thomas Klemm Wrote: Or if we start with two points whose function values have opposite signs and then evaluate the function at the midpoint. Not really, because in my idea no needed to bracket the root. The algorithm identifies the right direction, then make one step in that direction. When the sign is changed, the direction changes and the stepsize is reduced. Repeat until the stepsize is small enough (when the sign is changed). First I wrote in similar way an extremum finder for CASIO fx4000P, later reduced the length and coded on fx3650P (see the attachment: variables: A: stop condition, the minimal stepsize, B: function previous value, C: step direction: +1/1, D(>0!) actual stepsize, X: actual x, Y: actual y, M: subroutine return address, the function: 150×9.81×0.48/(cos(x)+0.48×sin(x)) ). And of course, later I found totally same  implemented on a HP9100: I have never learned computer science, maybe my methods are documented and invented many years before, and of course, it is possible to upgrade them (or available more effective), but my goal is find a workable simplest method which can be coded immediately. Csaba 

02252019, 11:00 PM
(This post was last modified: 02252019 11:18 PM by Thomas Klemm.)
Post: #25




RE: Small Solver Program
(02252019 08:39 PM)Csaba Tizedes Wrote: Not really, because in my idea no needed to bracket the root. There's nothing wrong with a digitbydigit algorithm. However after the first change of sign you already bracket the root. So let's assume that this interval has length \(I\). And we want to bring it down to a small value, say \(\epsilon\). Digitbydigit method With each digit we divide the interval by 10. After \(k\) digits we end up with: \(\epsilon=\frac{I}{10^k}\) Thus \(k=\frac{\log I  \log \epsilon}{\log 10}\). But on average we need 5 steps until we can decide on the digit. All in all we use about \(\frac{5}{\log 10}(\log I  \log \epsilon)\) steps. Binary search We divide the interval consecutively by 2. After \(k\) bisections we end up with: \(\epsilon=\frac{I}{2^k}\) Thus \(k=\frac{\log I  \log \epsilon}{\log 2}\). We only need 1 step to decide which interval to keep. All in all we use \(\frac{1}{\log 2}(\log I  \log \epsilon)\) steps. Conclusion Now compare \(\frac{5}{\log 10} \approx 2.1715\) with \(\frac{1}{\log 2} \approx 1.4427\) to see that the binary search is about \(1.5051\) times faster than the digitbydigit method. The only advantage to use a digitbydigit method is when calculating the function value of the next step can be done easier based on the function values of the previous steps. Therefore, it's used in calculating the square root or the trigonometric functions (CORDIC). However, I do not see any way to do this with an arbitrary function like the one given in the video: \(f(x)=\left (\frac{e \ln(\pi)}{\ln(x)} \right )^xe^\pi\) Cheers Thomas 

11032019, 03:14 PM
(This post was last modified: 11072019 01:07 AM by Albert Chan.)
Post: #26




RE: Small Solver Program
(02162019 03:58 AM)Thomas Klemm Wrote:(02152019 09:16 PM)Albert Chan Wrote: Amazingly, my rate formula is same as Aitken extrapolation formula ! Except for r=1, rate formula work even if iterations diverges (r > 1) Example, solve x = f(x) = 5  x^3, guess = 1.5 Above setup, iterations diverges: 1.5 → 1.625 → 0.708984375 → 4.643622734 ... First 3 numbers, r = (0.708984375  1.625) / (1.624  1.5) ≈ 7.328 Iterate with x = g(x) = p(5  x^3) + (1p)x, with p = 1/(1r) ≈ 12.0% : 1.5 1.515000000 1.515928095 1.515977481 1.515980083 1.515980220 1.515980227 1.515980228 ← converged We can speed up the convergence by using r = f'(guess) = 3 * 1.5^2 = 6.75 Iterate with same g(x), but p = 1/(1r) ≈ 12.9% : 1.5 1.516125000 1.515977551 1.515980277 1.515980227 1.515980228 ← converged 

11102019, 07:02 PM
Post: #27




RE: Small Solver Program
(11032019 03:14 PM)Albert Chan Wrote: Except for r=1, rate formula work even if iterations diverges (r > 1) Proof: Assumed rate r = Δx_{i+1} / Δx_{i} = constant, we extrapolate for the converged value. if r < 1: x_{∞} = x_{0} + (x_{1}  x_{0}) + (x_{2}  x_{1}) + ... = x_{0} (1 + r + r² + ...) = x_{0}/(1r) if r > 1: we can visualize it as iterate backintime, with rate of 1/r x_{0}  x_{∞} = (x_{0}  x_{1}) + (x_{1}  x_{2}) + (x_{2}  x_{3}) + ... = x_{0} (1 + 1/r + 1/r² + ...) x_{∞} = x_{0} (1  1/(11/r)) = x_{0}/(1r) x_{∞} and x_{∞} have the same expression ! x_{0}/(1r) = (x_{1}  (x_{1}  x_{0})) / (1r) = (f(x_{0})  r x_{0}) / (1r) rate iteration formula: x ← (f(x)  r x) / (1r) 

12012019, 12:13 AM
(This post was last modified: 01042020 09:58 PM by Albert Chan.)
Post: #28




RE: Small Solver Program
(11102019 07:02 PM)Albert Chan Wrote: rate iteration formula: x ← (f(x)  r x) / (1r) The formula actually approximate Newton's method. Let g(x) = x  f(x) → solving for x=f(x) is same as solving g(x)=0 g'(x) = 1  f'(x) ≈ 1  r x ← (f(x)  r x) / (1r) ≈ ((x  g(x))  (1  g'(x))*x) / g'(x) = x  g(x) / g'(x) 

01042020, 07:49 PM
Post: #29




RE: Small Solver Program
Just realized rate formula is equivalent to secant's method !
Post #6 showed rate formula is equivalent to Aitken's formula. Post #27 showed we can visualize sequence running backwards, getting the same formula. Assuming iteration formula f, with b = f(a), c = f(b): \(Aitken(a,b,c) = Aitken(c,b,a) = a  {(ab)^2 \over (ab)(bc)}\) If we make 2 points: \((x_1,y_1) = (a,ab) \;,\; (x_2, y_2) = (b,bc)\), we have: \(Aitken(a,b,c) = x_1  y_1 \left({x_1  x_2 \over y_1  y_2}\right)\) This is the formula for secant's method, interpolate for y = x  f(x) = 0 QED 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)