(11C) CUBIC EQUATION

[Dieter]

(06-11-2018 08:19 AM)Gamo Wrote:  After I took off the three line and test it somehow program keep running none stop.
So I put that three line back and it work.

Sure. These three lines are required, if you remove them the program will not work. The lines are missing in the 29C listing, that's an error.

(06-11-2018 08:19 AM)Gamo Wrote:  Seem like after adding the [CLx] [STO] 4 to [LBL] A somewhat the program run a bit faster than before.

This can be expected. But it still is not running very fast. Here is why:

I tried to understand how the program works. The general method is obvious: Determine one (real) root of the cubic equation, then calculate the other two from this. The first root is calculated by iteration. Here the program uses a simple but not very effective method: the value is determined digit by digit!

Here is an example: For the first root of x³–x²–x+0,5 at –0,8546376797 the program generates the following sequence:

Code:

-0,1
-0,2
-0,3
-0,4
-0,5
-0,6
-0,7
-0,8
-0,9
-0,89
-0,88
-0,87
-0,86
-0,85
-0,851
-0,852
-0,853
-0,854
-0,855
-0,8549
-0,8548
-0,8547
-0,8546
-0,85461
-0,85462
-0,85463
-0,85464
-0,854639
-0,854638
-0,854637
-0,8546371
-0,8546372
-0,8546373
-0,8546374
-0,8546375
-0,8546376
-0,8546377
-0,85463769
-0,85463768

The iteration exits if |f(x)| drops below 1E-4 (or 1E-8 in my version below) or if the last two approximations match (because a change beyond the 10th digit would occur). For the x=–0,8546... example this requires about 40 iterations for f(x)<1E–8.

On the other hand, if the root has only a few significant digits, for instance if it's an integer like 1 or 3 in the first example, the iteration finds the first root quite fast. But if the value needs 9 or 10 digits to get f(x) below 1E–8 the average iteration count is 40 loops or more. That's why it takes so long. And that's probably also why the original error threshold was as high as 1E-4. This gives a mid-precision approximation to the first root in about half the time required for 1E–8. But since the other two roots are calculated from the first one they are also not exact, but just more or less accurate approximations. In other words: at least in the original version the program mostly does not give exact or near exact results, but solutions that may only have 5 correct digits.

Here is a modified version of the program.

• The program now handles the general case ax³+bx²+cx+d
• Unlike the original version, the constant term now may also be zero.
• Input is a [ENTER] b [ENTER] c [ENTER] d
• The f(x) threshold for the first root has been reduced to 1E-8.
  This way the returned roots are more accurate.
• The program now returns the solutions directly.
  There is no need to press different keys depending on the sign of D.
• Because of this the discriminant D is no longer displayed.
• The 2nd/3rd complex roots are displayed as "imaginary part" (PAUSE), real part.
• The 2nd/3rd real roots are displayed with R/S.
  

The listing:

Code:

001  LBL A     // input: a ENTER b ENTER c ENTER d
002  STO 3
003  R↓
004  STO 2
005  R↓
006  STO 1
007  R↓
008  STO÷1
009  STO÷2
010  STO÷3
011  EEX
012  CHS
013  8         // lower error threshold than the original 1E-4
014  STO 0
015  CLx
016  STO 4
017  RCL 3
018  STO 6
019  ABS
020  x≠0?      // TEST 0 on 15C
021  STO÷6
022  RCL 1
023  ABS
024  +
025  RCL 0
026  LBL 0
027  1
028  0
029  x
030  x≤y?
031  GTO 0
032  STO 7
033  LBL 2
034  1
035  0
036  STO÷7
037  RCL 6
038  CHS
039  STO 6
040  LBL 1
041  RCL 7
042  RCL 6
043  x
044  RCL 4
045  +
046  STO 5
047  RCL 4
048  x=y?      // TEST 5 on 15C
049  GTO 3
050  RCL 5
051  RCL 1
052  +
053  RCL 5
054  x
055  RCL 2
056  +
057  RCL 5
058  x
059  RCL 3
060  +
061  STO 8
062  ABS
063  RCL 0
064  x>y?      // TEST 7 on 15C
065  GTO 3
066  RCL 5
067  STO 4
068  RCL 8
069  RCL 6
070  x
071  x<0?      // TEST 2 on 15C
072  GTO 1
073  GTO 2
074  LBL 3
075  RCL 5
076  R/S       // first real root
077  RCL 5
078  RCL 1
079  +
080  CHS
081  STO 1
082  2
083  ÷
084  ENTER
085  X^2
086  RCL 1
087  CHS
088  RCL 5
089  x
090  RCL 2
091  +
092  STO 3
093  -         // = discriminant D
094  x<0?      // TEST 2 on 15C
095  GTO 4     // for D<0 go to complex roots
096  √X
097  X<>Y
098  ABS
099  +
100  RCL 1
101  ENTER
102  ABS
103  ÷
104  x
105  STO 5
106  R/S       // second real root
107  RCL 3
108  RCL 5
109  ÷
110  RTN       // third real root
111  LBL 4
112  CHS       // 2nd & 3rd complex roots
113  √X
114  PSE       // show imaginary part
115  X<>Y
116  RTN       // show real part

Usage:

Enter the four coefficients of the equation ax³ + bx² + cx + d:
a [ENTER] b [ENTER] c [ENTER] d

Press [A]
This returns the first root, calculated by iteration.

Press [R/S]
- If the other roots are complex: "imaginary part" (PAUSE) and real part are returned.
- Otherwise the second real root is returned, [R/S] then displays the third one.

Examples:

x³ – x² – x + 1/2
1 [ENTER] 1 [CHS] [ENTER] 1 [CHS] [ENTER] 0,5
  [A]  => –0,854637680
[R/S] =>   1,451605963
[R/S] =>   0,403031717

x³ + x² – x
1 [ENTER] 1 [ENTER] 1 [CHS] [ENTER] 0
  [A]  =>   0,000000000
[R/S] => –1,618033989
[R/S] =>   0,618033989

2x³ + 4x² + 4x + 2
2 [ENTER] 4 [ENTER] 4 [ENTER] 2
  [A]  => –1,000000000
[R/S] => "0,866025404"  –0,500000000
So that's –0,5 ±0,866 i.

Dieter