Untar the package of C code on the class website named 'sysroots.tar'. This will create a subdirectory named 'sysroots' under the current directory. Change into that 'sysroots' subdirectory.
In this section we will use the multidimsional version of the Newton-Raphson algorithm to find the roots of a system of nonlinear equations:
f(x,y) = x² + y² - 1 = 0The locus of points in the x,y plane that solve each of these two functions is shown in the following plot, where the red curve is for f and the green curve is for g:
Note that there are two roots of the system that we should be able to converge onto.
Now in a terminal window, prepare to edit the C source code file by starting a gedit process running in the background with:
gedit nonlin_sys.c &
Note the ampersand character '&' appended to this command line. This will cause the command to be run in the background but return to a bash shell prompt to allow bash shell commands for compiling and running C code to continue to be run in the terminal. This command starts a gedit process to edit the C source code, and you should see a separate gedit window pop up with the nonlin_sys.c source code ready for editing. Read through this short program and see how the main program accepts two command line arguments and then calls the function 'sys_newton_raphson', which is supplied by the system root finding tool code 'sys_roots.c'. At what point will the iterative solution be initialized? How many maximum iterations are allowed? At what tolerance will the iterations be terminated?
The function call to sys_newton_raphson passes two pointers to functions that compute the system of functions to be solved and the corresponding Jacobian matrix. Find those functions at the top of the nonlin_sys.c source file and complete the lines indicated with the comments so that the desired system and its Jacobian will be computed. Note that the sys_newton_raphson function makes use of the gauss_pivotmax() function in the lineq.c tool code to converge on the roots of the system.
Save your completed source code file and compile it linked together with both the sys_roots.c and the lineq.c code with the shell command:
gcc -Wall nonlin_sys.c sys_roots.c lineq.c -o nonlin_sys -lm
Now at the shell prompt, run the compiled program and note its help output line:
./nonlin_sys
This is a reminder that the program expects your starting guesses for an x value and a y value to initialize the algorithm. Since we have seen on the plot above that one root can be expected somewhere in the vicinity of x=1 and y=1, rerun the program giving it that starting guess. How many iterations are required to converge to the root for the given tolerance? Rerun the program with another starting guess designed to let it converge to the other system root. What happens if you supply a starting guess of x=0 and y=0? Why?
It turns out that this system of equations is a simple algebraic relationship that can be easily solved by the quadratic formula. Solve it analytically and compare your answers to the values the Newton-Raphson algorithm converged on.
Now consider the case when the analytical derivation of the partial derivatives of the system functions is difficult, and approximating them with numerical derivatives is a good choice. Perhaps we have this system of equations to solve:
f(x,y) = 4tanh(x) - y = 0The locus of points in the x,y plane that solve each of these two functions is shown in the following plot, where the red curve is for f and the green curve is for g:
We should be able to converge on the single root of this system. Now back in a terminal window, prepare to edit the C source code file that contains the numerical approximation to the Jacobian matrix calculation by starting a gedit process running in the background with:
gedit nonlin_sys2.c &
Complete the source code lines needed to compute the values of the two functions to be solved. Save your completed source code file and compile it linked together with both the sys_roots.c and the lineq.c code with the shell command:
gcc -Wall nonlin_sys2.c sys_roots.c lineq.c -o nonlin_sys2 -lm
Now at the shell prompt, run the compiled program and note its help output line:
./nonlin_sys2
This is a reminder that the program expects your starting guesses for an x value and a y value to initialize the algorithm. Since we have seen on the plot above that one root can be expected somewhere in the vicinity of x=1 and y=2, rerun the program giving it that starting guess. Check the converged values with the plot above.