Sunday, February 12, 2017

Computer Solutions/Integer Solutions of Nonlinear Systems of Equations


Jsun Yui Wong

The two computer programs listed below seek to solve a continuous case and an integer case of two nonlinear systems of equations, respectively.  The first case comes from page 44 of Greenspan and Casulli [4] and is as follows:  

         -10 * X(1) + 5 * X(2) - EXP(X(1)),

         5 * X(1) - 10 * X(2) + 5 * X(3) - EXP(X(2)),

         5 * X(2) - 10 * X(3) + 5 * X(4) - EXP(X(3)),

         5*X(3) - 10 * X(4) - EXP(X(4))=0.


0 DEFDBL A-Z

3 DEFINT J, K

4 DIM X(342), A(342), L(333), K(333)


12 FOR JJJJ = -32000 TO 32000

    14 RANDOMIZE JJJJ

    16 M = -1D+317

    91 FOR KK = 1 TO 4


        94 A(KK) = -10 + FIX(RND * 21)



    95 NEXT KK
    128 FOR I = 1 TO 500000



        129 FOR K = 1 TO 4


            131 X(K) = A(K)
        132 NEXT K
        155 FOR IPP = 1 TO FIX(1 + RND * 3)

            181 B = 1 + FIX(RND * 4)


            182 REM IF RND < -.1 THEN 183 ELSE GOTO 189

            183 R = (1 - RND * 2) * A(B)
            186 X(B) = A(B) + (RND ^ 3) * R

            188 GOTO 191
            189 IF RND < .5 THEN X(B) = A(B) - FIX(1 + RND * 1.99) ELSE X(B) = A(B) + FIX(1 + RND * 1.99)


        191 NEXT IPP



        291 FOR j44 = 1 TO 4
            293 IF X(j44) > 30 THEN 1670



        299 NEXT j44



        333 X(3) = (10 * X(4) + EXP(X(4))) / 5

        334 IF X(3) > 30 THEN 1670


        336 N96 = 5 * X(1) - 10 * X(2) + 5 * X(3) - EXP(X(2))
        338 N97 = 5 * X(2) - 10 * X(3) + 5 * X(4) - EXP(X(3))


        344 N98 = -10 * X(1) + 5 * X(2) - EXP(X(1))

        1335 P = -ABS(N96) - ABS(N97) - ABS(N98)



        1499 IF P <= M THEN 1670
        1657 FOR KEW = 1 TO 4


            1658 A(KEW) = X(KEW)
        1659 NEXT KEW
        1661 M = P

    1670 NEXT I
    1888 IF M < -.000001 THEN 1999

    1917 PRINT A(1), A(2), A(3), A(4), M, JJJJ

1999 NEXT JJJJ

This computer program was run with qb64v1000-win [8]. Copied by hand from the screen, the computer program’s complete output through JJJJ=-31982 is shown below:

-.282813848003122        -.4148956819788955      -.4148956819702162
-.2828138479984057      -3.60628082951564D-10      -31985

-.2828138480373464        -.414895681997618      -.4148956820101716
-.2828138480169833        -7.031838999971285D-10      -31983

-.2828138480379652        -.4148956820071726      -.4148956819847773
-.282813848005176         -6.283482067992452D-10      -31982

Above there is no rounding by hand; it is just straight copying by hand from the screen.

On a personal computer with a Pentium Dual-Core CPU E5200 @2.50GHz, 2.50 GHz, 960 MB of RAM and with qb64v1000-win [8], the wall-clock time through JJJJ=-31982 was 15 seconds, not including "Creating .EXE file..." time.


The second case--shown below--comes from page 47 of Greenspan and Casulli [4]:

        X(1) +X(2) + X(3) + X(4) =0,

        LOG(X(1)) + X(2) + LOG(X(3)) + X(4) = -2,

        X(1) * X(2) * X(3) - X(4) =0,

        X(1) - EXP(X(2)) + X(3) + EXP(X(4)) =2.
 

0 DEFDBL A-Z

3 DEFINT J, K

4 DIM X(342), A(342), L(333), K(333)


12 FOR JJJJ = -32000 TO 32000


    14 RANDOMIZE JJJJ

    16 M = -1D+317

    91 FOR KK = 1 TO 4

        94 A(KK) = -50 + FIX(RND * 101)


    95 NEXT KK
    128 FOR I = 1 TO 1000



        129 FOR K = 1 TO 4


            131 X(K) = A(K)
        132 NEXT K
        155 FOR IPP = 1 TO FIX(1 + RND * 3)

            181 B = 1 + FIX(RND * 4)


            182 IF RND < -.1 THEN 183 ELSE GOTO 189

            183 R = (1 - RND * 2) * A(B)
            186 X(B) = A(B) + (RND ^ 3) * R

            188 GOTO 191
            189 IF RND < .5 THEN X(B) = A(B) - FIX(1 + RND * 1.99) ELSE X(B) = A(B) + FIX(1 + RND * 1.99)


        191 NEXT IPP



        1001 X(1) = -X(2) - X(3) - X(4)

        1002 IF X(1) < .01 THEN 1670
        1003 IF X(3) < .01 THEN 1670



        1004 N96 = LOG(X(1)) + X(2) + LOG(X(3)) + X(4) + 2

        1007 N97 = X(1) * X(2) * X(3) - X(4)

        1009 N98 = X(1) - EXP(X(2)) + X(3) + EXP(X(4)) - 2


        1335 P = -ABS(N96) - ABS(N97) - ABS(N98)



        1499 IF P <= M THEN 1670
        1657 FOR KEW = 1 TO 4


            1658 A(KEW) = X(KEW)
        1659 NEXT KEW
        1661 M = P

    1670 NEXT I
    1888 IF M < -1 THEN 1999

    1917 PRINT A(1), A(2), A(3), A(4), M, JJJJ

1999 NEXT JJJJ

This computer program was run with qb64v1000-win [8]. Copied by hand from the screen, the computer program’s complete output through JJJJ=-31911 is shown below:

1      -1      1      -1      0
-31982

1      -1      1      -1      0
-31947

1      -1      1      -1      0
-31936

1      -1      1      -1      0
-31926

1      -1      1      -1      0
-31923

1      -1      1      -1      0
-31917

1      -1      1      -1      0
-31911

Above there is no rounding by hand; it is just straight copying by hand from the screen.

On a personal computer with a Pentium Dual-Core CPU E5200 @2.50GHz, 2.50 GHz, 960 MB of RAM and with qb64v1000-win [8], the wall-clock time through JJJJ=-31911 was 2 seconds, not including "Creating .EXE file..." time.

Acknowledgment


I would like to acknowledge the encouragement of Roberta Clark and Tom Clark.

References
[1] R. Burden, J. Faires, A. Burden, Numerical Analysis, Tenth Edition. Cengage Learning, 2016.
[2] R. Burden, J. Faires, Numerical Analysis, Sixth Edition. Brooks/Cole Publishing Company, 1996.
[3] R. Burden, J. Faires, Numerical Analysis, Third Edition. PWS Publishers, 1985.
[4] D. Greenspan, V. Casulli, Numerical Analysis for Applied Mathematics, Science, and Engineering. Addison-Wesley Publishing Company, 1988
[5] L. W. Johnson, R. D. Riess, Numerical Analysis, Second Edition. Addison-Wesley Publishing Company, 1982
[6] Microsoft Corp. BASIC, second edition (May 1982), Version 1.10. Boca Raton, Florida: IBM Corp., Personal Computer, P. O. Box 1328-C, Boca Raton, Florida 33432, 1981.
[7] William H. Mills, A System of Quadratic Diophantine Equations, Pacific Journal of Mathematics, 3 (1953), pp. 209-220.
[8] Wikipedia, QB64, https://en.wikipedia.org/wiki/QB64.

No comments:

Post a Comment