Proof.
For a complex calculation, we only need to calculate its Re and Imaccurately. Let . In order toprove that there exist a
1)For addition , our algorithm is
Now let's analyse its error.
For the relative accuracy, let
Re and Im have high relative accuracy.
2)For substraction , our algorithm is
Now let's analyse its error
For the relative accuracy, let
Re and Im have high relative accuracy.
3)For multiplication , our algorithm is
Now let's analyse its error
For relative accuracy
When
4)For division
Now let's analyse its error, we only analyse thecase
For relative accarucy,
Cause we won't multiple two big number, the problemof overflow won't happend.
This is properly tested on python code as the following.
>>> def complex_divide(a1, a2, b1, b2):Question. Prove Lemma 1.3. Let
Proof.
We only prove the case in
The matrix form of inner product. Let
We shall show that A is s.p.d.,
which means
If A is s.p.d, we shall show that
1)
2).
3)
4)for A is h.p.d,
Question 1.14Question. Prove lemma 1.5:
Proof.
Question 1.15Question. Prove Lemma 1.6:
An operator norm is a matrix norm.
Proof.
1)
If
2)
3)
Question 1.16Question. Prove all parts except 7 of Lemma 1.7. Hint for part 8:Use the fact that if
Part 1. Prove that
Proof. For operator norm, this can be soon obtained by
if for any nonzero vector
For matrix Frobenius norm, if we rewrite
The last step is to sum the square of all the inequalities above,i.e.
Part 2. Prove that
Proof. For operator norm,
For Frobenius norm, this can be soon obtained byrewrite
Part 3. Prove that the max norm and Frobenius norm are not operatornorm.
Proof. For max norm, we may show that it violates part 2 ofLemma 1.7 by the counter example,
It canbe easily shown thatFrobenius norm is not operator norm, as
Part 4.
Proof. Only to Prove when orthogonal, since unitary is the same.For operator norm induced by
For Frobenius norm, it is the Pythagorean theorem. if
Part 5. Prove that,the maximum absolute row sum.
Proof. this can be obtained by take
Part 6. Prove that,the maximum absolute column sum.
Proof. This is trivial from the view of functional analysis, as
Part 8. Prove that
Proof. This is the same as the proof of part 6.
Part 9. Prove that
Proof. Remember that
This is true for
Part 10. If
Proof. There exists
wherethe last inequality is a corollary from the inequality of means.There exists Part 11. If
Proof. This can be soon obtained from parts 6, 8, and 10.
Part 12. If
Proof. This can be soon obtained from parts 10 and 11.
Part 13. If
Proof. For the first part, find one, then define
The second part is true, since
Question 1.17Question. We mentioned that on a Cray machine the expression
Proof.
Now consider the case that we want to calculate
For any
As
Define Andobviously,
We want to prove that
...(Wait for completion. If you can do it, please contact us)
When
Question 1.18Question. Suppose that
if
Prove the following facts:
Barring overflow or underflow, the only round-off error committed inrunning the algorithm is computing
Thus, this program in effect simulates quadruple precision arithmetic,representing the true sum
Using this and similar tricks in a systematic way, it is possible toefficiently simulate all four basic floating-point operations inarbitrary precision arithmetic, using only the underlying floating pointinstructions and no "bit-fiddling". 128-bit arithmetic is implementedthis way on the IBM RS6000 and Cray (but much less efficiently on theCray, which does not have IEEE arithmetic).
Proof. We shall prove the following lemma: If floating point numbers
Proof of lemma. Remember we are in the binary system, and
Assume
When
When
When
For the case
Question 1.19Question. This question illustrates the challenges in engineeringhighly reliable numerical software. Your job is to write a program tocompute the two-norm
This algorithm is inadequate because it does not have the followingdesirable properties:
It must compute the answer accurately (i.e., nearly all the computeddigits must be correct) unless
It must be nearly as fast as the obvious program above in mostcases.
It must work on any "reasonable" machine, possibly including ones not running IEEE arithmetic. This means it may not cause an errorcondition unless
To illustrate the difficulties, note that the obvious algorithm failswhen
This routine is important enough that it has been standardized as aBasic Linear Algebra Subroutine, or BLAS, which should be availableon all machines.
Answer.
To avoid overflow, one can separate out the exponent part of a floating-point number, and compute the scaled squares. One result can be like
import numpy as npThis code is tested for
>>> norm2([2**-1024]) == 2**-1024Though for the convenience of coding, we cannot apply the built-infunction <frexp> in C code and implement it in a slow way, theperformance can be still evaluated. Assume the 'frexp' operation can beperformed in a proper way (in machine level, separate the exponent partand the tail part), and calculated as three operations, totally thisalgorithm cost
Further, the code is implemented in C code. (Please go see attachment ingitee)
https://gitee.com/j7168908jx/shuzhidaishu/blob/master/code/c1-19/code.c
Test results are given in table 11. norm2simple is the baseline algorithm given inthe question, norm2lapack is the one from LAPACK in Fortran code.norm2 and norm22 are the algorithms designed in the answer, wherenorm22 uses more efficient but probably harder to read operations, andmight not work with denormalized floating-point numbers. The vector
During the completion of the C code, conclusions are obtained. One ofthe most important for the problem is that sometimes it is not soexpensive to do the division, rather than using other techniques thatmight add more complexity to the code and the running time. Limited CPUstructure(register cache, variable storage) is to blame.
Table 1. Test result for question 1.19Iterations, length(Question. We will use a Matlab program to illustrate how sensitivethe roots of a polynomial can be to small perturbations in thecoefficients. Polyplot takes an input polynomial specified by its rootsr and then adds random perturbations to the polynomial coefficients,computes the perturbed roots, and plots them. The inputs are
The first part of your assignment is to run this program for the following inputs. In all cases choose m high enough that you get afairly dense plot but don't have to wait too long. m = a few hundred or perhaps 1000 is enough. You may want to change the axes of the plot if the graph is too small or too large.
Also try your own example with complex conjugate roots. Which rootsare most sensitive?
r=(1:10);e=1e-3,1e-4,1e-5,1e-6,1e-7,1e-8,
r=(1:20);e=1e-9,1e-11,1e-13,1e-15,
r=[2,4,8,16,...,1024];e=1e-1,1e-2,1e-3,1e-4.
The second part of your assignment is to modify the program tocompute the condition number c(i) for each root. In other words, arelative perturbation of e in each coefficient should change rootr(i) by at most about e*c(i). Modify the program to plot circlescentered at r(i) with radii e*c(i), and confirm that these circlesenclose the perturbed roots (at least when e is small enough thatthe linearization used to derive the condition number is accurate).You should turn in a few plots with circles and perturbed roots, andsome explanation of what you observe.
In the last part, notice that your formula for c(i) "blows up" if
Answer.
1.Here we let m =1000.
1)The case r=1:10, let e=1e-3,1e-4,1e-5,1e-6,1e-7,1e-8. See figures below.
Figure. When r=1:102)The case r=1:20,let e=1e-9,1e-11,1e-13,1e-15, See figures below.
Figure. When r=2:103)The case r=[2,4,8,...,1024],let e=1e-1,1e-2,1e-3,1e-4, See figures below.
Figure. When r=[2,4,8,...,1024]2.Let's start from the theoretical analysis. Let,
where
So
For our numerical experiment, we assume that
r=1:10,e=1e-7We can see that the circle in the middle is muchlarger than the circle on the edeg, but the circle of the roots 5 or 6is not the largest . It is easy to see that
derivative3.The same trick as 2, notice that k-th derivative of p(x) will be 0 if
See the code appendix c1-20
https://gitee.com/j7168908jx/shuzhidaishu/tree/master/code/c1-20
Question 1.21Question. Apply Algorithm 1.1, Bisection, to find the roots of
Answer.
The bisection is implemented as
import numpy as npand the great difference in result can be inherited from figure listed below.
result for question 1.21Modified bisection search code is listed below.
def bisect(roots, a, b, tol):Respectively, the outcome of code, when given different initial value,is shown in the figure below.
result for question 1.21有什麼問題歡迎在下方留言哦~ <