﻿ Euler Math Toolbox - Examples

# Newton-Algorithm for two Variables

We solve

>function f([x,y]) &= [x^2+y^2-10,x+y-1]
2    2
[y  + x  - 10, y + x - 1]

Euler can plot the solutions of both equations.

>plot2d(&f(x,y)[1],level=0,r=5); ...

The intersections are the solutions.

To start the Newton algorithm we need the Jacobian.

>function Df([x,y]) &= jacobian(f(x,y),[x,y])
[ 2 x  2 y ]
[          ]
[  1    1  ]

Now, there is the newton2 function, which does the iteration.

It needs to function f(v) and Df(v). That is, why we define f and Df with f([x,y]) and Df([x,y]). Those functions can be used with vectors or two elements.

>newton2("f","Df",[-3,3])
[-1.67945,  2.67945]

We want to simulate the Newton iteration step by step. So we define the iterating function.

>function fiter ([x,y]) &= [x,y]-invert(Df(x,y)).f(x,y)
[    2    2                            ]
[ - y  - x  + 10   2 y (y + x - 1)     ]
[ -------------- + --------------- + x ]
[   2 x - 2 y         2 x - 2 y        ]
[                                      ]
[   2    2                             ]
[  y  + x  - 10   2 x (y + x - 1)      ]
[  ------------ - --------------- + y  ]
[   2 x - 2 y        2 x - 2 y         ]

>&factor(fiter(x,y))
[    2          2      ]
[ - y  + 2 y - x  - 10 ]
[ -------------------- ]
[      2 (y - x)       ]
[                      ]
[   2    2             ]
[  y  + x  - 2 x + 10  ]
[  ------------------  ]
[      2 (y - x)       ]

Maxima returns a column vector. We need a row vector as input and output. So we define an iteration function g.

>function g(v) := fiter(v)'

Then we can iterate using the iterate() function.

>iterate("g",[-3.2,3],10)
-3.2             3
-1.87419       2.87419
-1.68744       2.68744
-1.67946       2.67946
-1.67945       2.67945
-1.67945       2.67945
-1.67945       2.67945
-1.67945       2.67945
-1.67945       2.67945
-1.67945       2.67945
-1.67945       2.67945

Another starting point yields the other solution.

>iterate("g",[1,-2],10)
1            -2
3.16667      -2.16667
2.72396      -1.72396
2.67989      -1.67989
2.67945      -1.67945
2.67945      -1.67945
2.67945      -1.67945
2.67945      -1.67945
2.67945      -1.67945
2.67945      -1.67945
2.67945      -1.67945

Of course, we can also solve this example symbolically.

>&solve(f(x,y)), %()
1 - sqrt(19)      sqrt(19) + 1
[[y = ------------, x = ------------],
2                 2
sqrt(19) + 1      1 - sqrt(19)
[y = ------------, x = ------------]]
2                 2

[-1.67945,  2.67945,  2.67945,  -1.67945]

There is also the Broyden algorithm, which does not need the derivative. It is a generalization of the Secant algorithm for more then one variable.

>broyden("f",[-3,2])
[-1.67945,  2.67945]

With the interval Newton method, we even get a proved inclusion of the zero.

>inewton2("f","Df",[-3,1])
[ ~-1.679449471770339,-1.679449471770335~,
~2.679449471770335,2.679449471770338~ ]

Examples