::::::::::::::
as2.sol
::::::::::::::
1. See files as2q1.c, as2q1.out
For (a) and (b), the inaccurate results are due to GE with no pivoting
being an unstable algorithm; a large multiplier occurs in these 2 cases.
For (c) and (d), all 3 algorithms produce inaccurate results (should be
x_1 = 1, x_2 = 10^5 and 10^6, resp.) because of an ill-conditioned
matrix that is nearly singular and roundoff errors in representing
numbers in base 2; note that the ranking of the 3 algorithms in order
of closeness to true answer is different in these 2 cases.
2. (a) Let d_i be the diagonal elements of D_1 and e_i be the diagonal
elements of D_2. Then b_ij = d_i*e_j*a_ij. I.e., B is obtained
from A by first multiplying the ith row of A by d_i, i = 1,...,n,
and then multiplying the jth column of D_1*A by e_j, j = 1,...,n.
(b) Let C = AP where P is a permutation matrix. Then b = A*x =
C*P^(-1)*x = C*y where y = P^(-1)*x = P^T*x. I.e., y is obtained
from x by permuting the elements of x in the same way as the columns
of A are permuted to get C.
Alternative: Let a_j and c_j = a_p(j) be the columns of A and C,
where p(1),p(2),...,p(n) is a permutation of 1,2,...,n.
Then b = x_1*a_1 + x_2*a_2 + ... + x_n*a_n
= x_p(1)*a_p(1) + x_p(2)*a_p(2) + ... + x_p(n)*a_p(n)
= y_1*c_1 + y_2*c_2 + ... + y_n*c_n
where y_j = x_p(j), j = 1,...,n.
3. A = [ 2 -2 3 3 ] b = [ 31 ]
[ -2 6 -4 2 ] [ -74 ]
[ 4 -4 8 -8 ] [ 128 ]
[ 1 5 -1 1 ] [ -37 ]
(a) Step 1: Pivot row is r_1 = 3; modified matrix with multipliers is
[ 4 -4 8 -8 ]
[ -1/2 4 0 -2 ]
[ 1/2 0 -1 7 ]
[ 1/4 6 -3 3 ]
Step 2: Pivot row is r_2 = 4; modified matrix with multipliers is
[ 4 -4 8 -8 ]
[ 1/4 6 -3 3 ]
[ 1/2 0 -1 7 ]
[ -1/2 2/3 2 -4 ]
Step 3: Pivot row is r_3 = 4; modified matrix with multipliers is
[ 4 -4 8 -8 ]
[ 1/4 6 -3 3 ]
[ -1/2 2/3 2 -4 ]
[ 1/2 0 -1/2 5 ]
The P, L, U matrices in PA = LU decomposition are:
P = [ 0 0 1 0 ] L = [ 1 0 0 0 ] U = [ 4 -4 8 -8 ]
[ 0 0 0 1 ] [ 1/4 1 0 0 ] [ 0 6 -3 3 ]
[ 0 1 0 0 ] [ -1/2 2/3 1 0 ] [ 0 0 2 -4 ]
[ 1 0 0 0 ] [ 1/2 0 -1/2 1 ] [ 0 0 0 5 ]
The matrices such that U = M_3*P_3*M_2*P_2*M_1*P_1*A are:
P_1 = [ 0 0 1 0 ] P_2 = [ 1 0 0 0 ] P_3 = [ 1 0 0 0 ]
[ 0 1 0 0 ] [ 0 0 0 1 ] [ 0 1 0 0 ]
[ 1 0 0 0 ] [ 0 0 1 0 ] [ 0 0 0 1 ]
[ 0 0 0 1 ] [ 0 1 0 0 ] [ 0 0 1 0 ]
M_1 = [ 1 0 0 0 ] M_2 = [ 1 0 0 0 ] M_3 = [ 1 0 0 0 ]
[ 1/2 1 0 0 ] [ 0 1 0 0 ] [ 0 1 0 0 ]
[ -1/2 0 1 0 ] [ 0 0 1 0 ] [ 0 0 1 0 ]
[ -1/4 0 0 1 ] [ 0 -2/3 0 1 ] [ 0 0 1/2 1 ]
(b) Forward substitution applied to Lz = Pb = [ 128 -37 -74 31 ]^T
yields z = [ 128 -69 36 -15 ]^T.
Back substitution applied to Ux = z yields x = [ -2 -4 12 -3 ]^T.
Forward substitution applied to Lw = Px = [ 12 -3 -4 -2 ]^T
yields w = [ 12 -6 6 -5 ]^T.
Back substitution applied to Uy = w yields y = [ -1 0 1 -1 ]^T.
4. We can interchange row_i and row_j as follows
row_i := row_i + (-1)*row_j
row_j := row_j + row_i
row_i := row_i + (-1)*row_j
row_i := (-1)*row_i
::::::::::::::
as2q1.c
::::::::::::::
#include
main()
{
int k,nsys;
double a[2][2],b[2],xi[2],xn[2],xp[2];
void ainv2(),genp2(),gepp2();
scanf("%d", &nsys);
for (k = 1; k <= nsys; k++) {
scanf("%lf %lf %lf", &a[0][0],&a[0][1],&b[0]);
scanf("%lf %lf %lf", &a[1][0],&a[1][1],&b[1]);
ainv2(a,b,xi);
genp2(a,b,xn);
gepp2(a,b,xp);
printf("A = %21.14e %21.14e b = %21.14e\n", a[0][0],a[0][1],b[0]);
printf(" %21.14e %21.14e %21.14e\n", a[1][0],a[1][1],b[1]);
printf("xi= %21.14e xn= %21.14e xp= %21.14e\n", xi[0],xn[0],xp[0]);
printf(" %21.14e %21.14e %21.14e\n", xi[1],xn[1],xp[1]);
printf("\n");
}
} /* main */
void ainv2(double a[2][2], double b[2], double x[2])
/* Given 2 by 2 matrix A, compute x = A^(-1)*b.
Input parameters: a, b. Output parameter: x. */
{
double det;
det = a[0][0]*a[1][1] - a[1][0]*a[0][1];
x[0] = (a[1][1]*b[0] - a[0][1]*b[1])/det;
x[1] = (a[0][0]*b[1] - a[1][0]*b[0])/det;
} /* ainv2 */
void genp2(double a[2][2], double b[2], double x[2])
/* Use Gaussian elimination with no pivoting
to solve 2 by 2 linear system A*x = b.
Input parameters: a, b. Output parameter: x. */
{
double a11,b1,m;
m = a[1][0]/a[0][0];
a11 = a[1][1] - m*a[0][1];
b1 = b[1] - m*b[0];
x[1] = b1/a11;
x[0] = (b[0] - a[0][1]*x[1])/a[0][0];
} /* genp2 */
void gepp2(double a[2][2], double b[2], double x[2])
/* Use Gaussian elimination with partial pivoting
to solve 2 by 2 linear system A*x = b.
Input parameters: a, b. Output parameter: x. */
{
double a11,b1,m;
int r0,r1;
if (fabs(a[0][0]) >= fabs(a[1][0]))
{ r0 = 0; r1 = 1; }
else { r0 = 1; r1 = 0; }
m = a[r1][0]/a[r0][0];
a11 = a[r1][1] - m*a[r0][1];
b1 = b[r1] - m*b[r0];
x[1] = b1/a11;
x[0] = (b[r0] - a[r0][1]*x[1])/a[r0][0];
} /* gepp2 */
::::::::::::::
as2q1.in
::::::::::::::
4
7.0e-6 1.23 1230.000007
7.89e6 4.56 7.89456e6
7.0e-12 1.23 1230.000000000007
7.89e12 4.56 7.89000000456e12
2.0 3.0 300002.0
2.0000001 3.0 300002.0000001
2.0 3.0 3000002.0
2.0000001 3.0 3000002.0000001
::::::::::::::
as2q1.out
::::::::::::::
A = 7.00000000000000e-06 1.23000000000000e+00 b = 1.23000000700000e+03
7.89000000000000e+06 4.56000000000000e+00 7.89456000000000e+06
xi= 1.00000000000000e+00 xn= 9.99999981234266e-01 xp= 1.00000000000000e+00
1.00000000000000e+03 1.00000000000000e+03 1.00000000000000e+03
A = 7.00000000000000e-12 1.23000000000000e+00 b = 1.23000000000001e+03
7.89000000000000e+12 4.56000000000000e+00 7.89000000456000e+12
xi= 1.00000000000000e+00 xn= 1.03942251631192e+00 xp= 1.00000000000000e+00
1.00000000000000e+03 1.00000000000000e+03 1.00000000000000e+03
A = 2.00000000000000e+00 3.00000000000000e+00 b = 3.00002000000000e+05
2.00000010000000e+00 3.00000000000000e+00 3.00002000000100e+05
xi= 1.00000761464965e+00 xn= 1.00041893598973e+00 xp= 9.99992326135580e-01
9.99999997207093e+04 9.99999997207093e+04 1.00000000005116e+05
A = 2.00000000000000e+00 3.00000000000000e+00 b = 3.00000200000000e+06
2.00000010000000e+00 3.00000000000000e+00 3.00000200000010e+06
xi= 9.99619563576836e-01 xn= 1.00535351317376e+00 xp= 9.99923261410584e-01
9.99999996430991e+05 9.99999996430991e+05 1.00000000005116e+06