Description

Ian Nathaniel J. Lapuz ijlapuz@up.edu.ph
Elaine T. Pajarillo etpajarillo@up.edu.ph
Junius Wilhelm G. Bueno
Gov. Pack Road, Baguio City, Baguio 2600
November 3, 2023
1. Consider the following pentadiagonal matrix,
a. Make a Python function block that prints matrix A for different n.

In this code, we defined the ‘block(n)’ function to create the Matrix A based on the input value n and the given initial values for a, b, c, d, and e. First, we initialize an n-by-n matrix A with elements set to zero, and inputted the given initial values for variables a, b, c, d, and e from the problem. We created the matrix A by assigning the values of c[i] as the diagonal elements. We apply the conditions to populate the nondiagonal elements: if i>0 or if i is not on the first row or column, we assign the values of b[i] as the upper diagonal elements of the matrix and d[i-1] as the elements on the lower diagonal. Furthermore, if i>1 or i is at least two rows or columns away from the diagonal, we assign a[i] to the elements in the second upper diagonal, and e[i] to the elements in the second lower diagonal. After completing the process, the code will return the Matrix A for any values of n.
b. Include in linearsys.py the column-wise Forward and Backward substitution methods and modify LUSolve to use the column-wise substitution method

This code contains the Columnwise Forward and Backward substitution methods, this code harbors much resemblance with its row-wise counterpart, with only the difference in the columnwise method iterating through columns first, and then iterating through the rows for each of the columns in order to compute for the solution. The Columnwise Forward Substitution Method will update the elements in order of columns, from the first to the last, while the Columnwise Backward Substitution Method will update the elements in the reverse order, from the last to the first. As for the modification of the LUSolve, we just called the column-wise substitution method instead of the rowise substitution method.
c. LU factorization of matrix A and x˜, y˜ ∈ ℝ that satisfies: Ax˜ = c, for n = 150
Solving for the Ax˜ = c, we initially stored the values of c given by the problem in an array of range 150. The resulting array of values saved in CVal to be used in computing for the solution. Then we used LUSolve to compute the solution by utilizing the previously computed A_matrix and CVal to be saved in solution1. Additionally, we also computed the relative error incurred by the solution by getting the norm of the CVal subtracted by the dot product of matrix A and the solution1 divided by the norm of the CVal. Hence, we have arrived at the output above for the computation of the solution x˜ and its corresponding relative error.
d. LU factorization of matrix A and x˜, y˜ ∈ ℝ that satisfies: Ay˜ = d, for n = 100
In solving for the Ay˜ = d, the same process above was performed. First saving the values of d given by the problem in an array of range 100 with the resulting array of values saved in DVal to be used in computing for the solution. Then LUSolve is used to compute for the solution using the previously computed A_matrix and DVal to be saved in solution2. We also computed the relative error that the solution incurred by getting the norm of the DVal subtracted by the dot product of matrix A and the solution2 divided by the norm of the DVal. Hence, we have arrived at the output above for the computation of the solution y˜ and the relative error.
2. Use the Jacobi method to find the solution of the linear system with initial iterate x0 = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], tolerance 10−15, and maximum iterate of 1000.
In order to use the Jacobi Method to find the solution of the given linear system, we have to first ensure that the main diagonal of our system has only non-zero elements. Putting the given system inside Matrix A with its corresponding vector b our initial inputs are given by:

From this, we can see that there are elements equal to zero in the main diagonal which prevents the implementation of the Jacobi Method, to resolve this, we need to perform partial pivoting on Matrix A to ensure a main diagonal with all non-zero elements. Below is the code that we used in solving for the solution of the linear system with the additional step of Partial Pivoting the Matrix.

The code starts with performing partial pivoting, which is used to ensure that the main diagonal of the matrix does not contain zero elements. In this process, a simple for loop checks if the pivot element is greater than the absolute value of the diagonal element at position [i][j]. If it is, rows are swapped to rearrange the matrix until the absolute value of the elements on the diagonal is greater than the absolute value of all the other elements. During the process of this row swapping, the corresponding values of the matrix A in vector b are also rearranged accordingly. Consequently, after completing the partial pivoting, both Matrix A and vector b are rearranged so that the main diagonal of the matrix does not contain any zero elements.
Following the partial pivoting of Matrix A, we output the newly rearranged matrix as ‘Sorted Matrix.’ Then, we proceed to apply the Jacobi Method to solve the linear system. This method can now be used because the initial condition of having all non-zero diagonal elements has been met. We call a function from the ‘linearsys.py’ file to compute the solution, the number of iterations, and the relative error and just output the results accordingly.
Upon running “Item2Imp.py”, the following results are yielded:
Solution:
x1 = -0.7823189693469743
x2 = 4.064637938693944
x3 = -4.0204353620613675
x4 = 1.52210128831629
x5 = -0.05442025766326386
x6 = 1.3698911594846728
x7 = 0.23300755219902242
x8 = 2.10242669924478
x9 = -0.08834962239004884
x10 = 1.7279125944024878
Number of Iterations: 243
Relative Error: 9.53223792001001e-16
3. Find the solution [a ∗ , b∗ , c∗ , d∗ ] to the nonlinear system of equations:
The following Python code “Item3Imp.py” is used to find the solutions of the given nonlinear systems of equations. It imports “linearsys.py” so that we can use the LUSolve() function.

This first part of the code contains the function newton(f, Jf, x, tol, maxit), the Python implementation of the Newton Method. It takes in the parameters f, Jf, tolerance, and max iterations. Also, for both items, we let x[0] = a, x[1] = b, x[2] = c, and x[3] = d, respectively, and convert the given functions to the said form for efficiency. with initial iterate x = [1, 1, 1, 1], tolerance 10
For item A, note that:
0 = 20[0][2] − 8[0][1] + 16[2]3 − 4[1][3] − 39 = 0
2
Solvingfor1 ===their 124−[partial03[0]2][0++]2derivates,[63[][11−]] [22+][ 1−2we]210[2+get:]7[−2][22+][3[23]][−0−]1116[3=]2=+ 0 07 = 0
3
/ ∂x[0] ∂x[1] ∂x[2] ∂x[3]
∂0 20[2] − 8[1] − 8[0] − 4[3] 20[0] + 48[2]2 − 4[1]
∂1 12 6 2 − 2
∂2 8[0] + 2[3]2 2[2] 2[1] − 10 4[0][3]

∂3 − 3[3] − 4[1] 7[3] − 3[0] + 7[2]
Therefore, we have:
00 = 20[2] − 8[1] 20 = 8[0] + 2[3]2
01 == 20−8[0][0+] − 4[32] = 2[2]
02 48[2]
03 = − 4[1] 23 = 4[0][3]
10 = 12 30 = − 3[3] 11 = 6 31 = − 4[1]
12 = 2 32 = 7[3]
We, then,13 hardcode= − 2 them into “Item3Imp.py”33under= −the3function[0] + 7af(x)[2] and aJf(x) so that
arrays of the inputted functions are returned. There is function aA(x) contains the same functions but without the b values, which will be used for the computation of the relative error.

To find the solution of Item A, we pass af(), aJf(), the initial iterate x = [1,1,1,1], the tolerance of 1e-14, and the max iteration of 100 to the newton() function. For the solving of the relative error, we make use of the imported np.linalg.norm(), aA(x), and the separated given solution “gA”.

Upon running “Item3Imp.py”, the following results are yielded for item A:
Item A:
Solution:
a* = 0.2669201878036756
b* = 1.5623911642242143
c* = 1.4480821232798882
d* = 2.2367767427745844
Error: 8.299304550897164e-17
Iterations: 7
[1, 1, 1], tolerance 10
For item B, note that:
Solvingfor01 ===their 3[−0[]partial[020]]− [1−]81+ derivates,(20[(1][[21+]]0−[.2we1]3))−2310get:−+π0=. 5 ( 0=[2 ]0) + 1. 06 = 0
2
/ ∂x[0] ∂x[1] ∂x[2]
∂0 3 [1]([1][2]) 2[0]
∂1 2[0] − 162([1] + 0. 1 ([2])
∂ −[0][1] −[0][1]
2 − [1] − [0] 20
Therefore, we have:
000102 === 3 [[21]](([[11]][[22]])) 122021 === −− ([[[012]]])−−[[00]][[11]]
10 = 2[0] 22 = 20
We, again,11 =hardcode − 162(them[1] +into0. 1“Item3Imp.py”) under the function bf(x) and bJf(x) so that
arrays of the inputted functions are returned. There is function bA(x) contains the same functions but without the b values, which will be used for the computation of the relative error.

For the solution of Item B, we pass bf(), bJf(), the initial iterate x = [1,1,1], the tolerance of 1e-14, and the max iteration of 100 to the newton() function. For the solving of the relative error, we make use of the imported np.linalg.norm(), bA(x), and the separated given solution
“gB”.
Python print(“Item B:”)
itemb=newton(bf,bJf,[1,1,1],1e-14,100)
gB=ls.np.array([0.5,-1.06,(3-(10*ls.np.pi))/3],float)#given solution array print(” Solution:”) forninrange(3):
print(f” {unk[n%4]}={itemb[0][n]}”)
gB = ls.np.array([0.5, -1.06, (3 – (10 * ls.np.pi))/3], float)
print(f”Error:{(ls.np.linalg.norm(bA(itemb[0])-gB))/ls.np.linalg.norm(gB)}”) print(f”Iterations:{itemb[2]} “)
Upon running “Item3Imp.py”, the following results are yielded for item B:
Item B:
Solution:
a* = 0.49999999999999994
b* = -1.657528818588421e-17
c* = -0.5235987755982988
Error: 1.8900455247697116e-16
Iterations: 8