# Shiyou's numerical analysis assignment

2020-11-08 08:04:34

### Numerical analysis assignment

It suddenly occurred to me that I could do numerical analysis , So I took my roommate's numerical analysis homework to practice , Write a blog to share . In fact, I'm a rookie who has complicated the program , There are many things that can be simplified , Because this novice has not finished other homework , I don't want to simplify , You can correct it by yourself .

# Numerical analysis on computer

First of all, let me talk about my doubts , For the first question python How to achieve to ln(x) Call the derivative directly ？ Is it a direct derivative of a polynomial after Taylor expansion ？ The second question is to use Newton Iterative method , It is required to find out the range of the initial value of iteration in which the convergent solution can be obtained , If we use iterative program to implement it, it will be more troublesome , Is there a better way to solve it ？ I hope there is a God who can help me to leave a message , Thank you very much .

# One 、 topic 1—— Approximate solution of logarithm

## 1. Title Description

** subject ：** There's a lot of laziness here , Just post the screenshot of the original question .

## 2.python Realization

Not much BB, The program is here ：

``````import numpy as np
from sympy import * # It is used in scientific calculation such as derivative integral
import math as m

x = Symbol('x')#x  Variable
t = Symbol('t')
y1 = 1/(1+x) # The formula
y2 = -1/(1+x) # The formula
y3 = 2/(1-x**2) # The formula

def func(m):
res = m
for j in range(1,m):
res *= j
return res

def ln_Tyalor(y):
Tl_expr = y * (t-x)
for i in range(1, 10):
fac = func(i+1)
f_n = diff(y, x, i)
Tl_expr += (f_n / fac)*(t-x)**(i+1)
return Tl_expr.subs({

x:0})

#print(ln_Tyalor(y1))
sexpr1 = ln_Tyalor(y1)
sexpr2 = ln_Tyalor(y2)
sexpr3 = ln_Tyalor(y3)
A = sexpr1.subs({

t:1}).evalf()
B = sexpr2.subs({

t:-1/2}).evalf()
C = sexpr3.subs({

t:1/3}).evalf()
print('ln2 Value ：', m.log(2, m.e))
print(' equation ln(1+x) Of 10 The result of Taylor expansion is ：', A,'\n',' The estimated error is ：', abs(m.log(2, m.e)-A))
print(' equation ln(1/(1+x)) Of 10 The result of Taylor expansion is ：', B,'\n',' The estimated error is ：', abs(m.log(2, m.e)-B))
print(' equation ln((1+x)/(1-x)) Of 10 The result of Taylor expansion is ：', C,'\n',' The estimated error is ：', abs(m.log(2, m.e)-C))
``````

## 3. Output results

``````ln2 Value ： 0.6931471805599453
equation ln(1+x) Of 10 The result of Taylor expansion is ： 0.645634920634921
The estimated error is ： 0.0475122599250246
equation ln(1/(1+x)) Of 10 The result of Taylor expansion is ： 0.693064856150793
The estimated error is ： 8.23244091517905e-5
equation ln((1+x)/(1-x)) Of 10 The result of Taylor expansion is ： 0.693146047390827
The estimated error is ： 1.13316911820593e-6

``````

# One 、 topic 2—— The root of an equation Newton Law

## 1. Title Description

** subject ：** Or the screenshots are convenient .

## 2.python Realization

This program is not well written , Because after writing ,exp（） Functions always report errors ： Integer data is not exp Object properties , I didn't realize my idea after I changed it , So be it , There's no time. , Let's do it ourselves .....

``````import numpy as np
from sympy import * # It is used in scientific calculation such as derivative integral
import math as m

def f(x):
return 3*x**2 - m.exp(x) # The equation has 3 A root , A preliminary estimate is that [-1,0],[0,1],[3,4]

def fdiff(x):
return 6*x - m.exp(x)

def g(x):
return 6*x - m.exp(x)# The equation has 3 A root , A preliminary estimate is that [0,1],[2,3]

def gdiff(x):
return 6 - m.exp(x)

a = float(input(' Please enter the lower bound of the calculation interval a( floating-point ）: '))
b = float(input(' Please enter the upper bound of the calculation interval b( floating-point ）: '))
c = float(input(' Please enter the iteration initial value ( floating-point ）: '))
n = input(' Please enter the function to be solved ,f representative f(x),g representative g(x): ')

if n =='f':

if f(a) * f(b)> 0:
print(' In this interval, the function has more than one root or has no root ')
elif f(a) * f(b) == 0:
print('f(a) = ', '%f'%f(a))
print('f(b) = ', '%f'%f(b))
else:
fcount = 0
y = c - f(c) / fdiff(c)
while (abs(c - y) >= 0.5e-9) & (fdiff(c) != 0):
x2 = c - f(c) / fdiff(c)
y = c
c = x2
fcount += 1
print(' The root interval given by the function is ：', [a, b])
print(" function f(x） Of Newton The number of iterations ：%f, function f(x） The root of the iterative calculation is :%.8f"%(fcount,c))

elif n =='g':

if g(a) * g(b)> 0:
print(' In this interval, the function has more than one root or has no root ')
elif g(a) * g(b) == 0:
print('g(a) = ', '%f'%g(a))
print('g(b) = ', '%f'%g(b))
else:
gcount = 0
y = c - g(c) / gdiff(c)
while (abs(c - y) >= 0.5e-9) & (gdiff(c) != 0):
x2 = c - g(c) / gdiff(c)
y = c
c = x2
gcount += 1
print(' The root interval given by the function is ：', [a, b])
print(" function g(x） Of Newton The number of iterations ：%f, function g(x） The root of the iterative calculation is :%.8f"%(gcount,c))
``````

## 3. Output results

Here's an explanation , Input [a, b] You want to determine if the interval has roots ;c Is the initial value of the iteration ;f representative f（x),g representative g(x); Please input these parameters by yourself .

`````` Please enter the lower bound of the calculation interval a( floating-point ）: -1.0
Please enter the upper bound of the calculation interval b( floating-point ）: 4.0
Please enter the iteration initial value ( floating-point ）: -3.0
Please enter the function to be solved ,f representative f(x),g representative g(x): f
The root interval given by the function is ： [-1.0, 4.0]
function f(x） Of Newton The number of iterations ：7.000000, function f(x） The root of the iterative calculation is :-0.45896227
``````

# One 、 topic 3——Hilbert The condition number of a matrix

## 1. Title Description

** subject ：** You'll see .

## 2.python Realization

``````import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['simhei']

n = int(input(' Please enter Hilbert The maximum dimension of a square matrix :' ))
def max_sum_rows(X):
sum_row_list1 = []
for i in range(len(X)):
count = 0
for j in range(len(X)):
count += abs(X[i][j])
sum_row_list1.append(count)
return max(sum_row_list1)

inf_fanshu = []
Hilbert_cond = []
for i in range(1, n+1):
X = 1./(np.arange(1, i+1) + np.arange(0, i)[:, np.newaxis])
invX = np.linalg.inv(X)
a1 = max_sum_rows(invX)
a2 = max_sum_rows(X)
inf_fanshu.append(a2)
H_cond = a1 * a2
Hilbert_cond.append(H_cond)

# Calculation 10,20……100 Infinite norm condition number of
Hilbert_cond_test = []
for j in range(n):
if (j+1)%10 == 0:
Hilbert_cond_test.append(Hilbert_cond[j])
print(' Generating dimensions from 10,20……100 Of Hilbert The row norm condition number of a matrix ：\n', Hilbert_cond_test)

plt.title('Hilbert The relationship between matrix dimension and logarithm of condition number ')
plt.plot((list(range(1,n+1))), np.log(Hilbert_cond),c ='b', marker='*',label=' Fit the curve ')
plt.legend()
plt.xlabel(' Matrix dimensions n')
plt.xticks(np.arange(0, n+1, 5))
plt.yticks(np.arange(0, 55, 5))
plt.ylabel('log(cond)')
plt.show()

# solve Hilbert The solution of the equation and the corresponding infinite condition number
r_A_A_acc_list = []
r_B_list = []
r_cond = []
r_B_cond = []
for i in range(1,n+1):
A = np.ones((i,1))*1
X = 1. / (np.arange(1, i + 1) + np.arange(0, i)[:, np.newaxis])
B = X@A
A_acc = np.linalg.inv(X)@B
r_A_A_acc = A - A_acc
r_B = B - X @ A_acc
r_A_A_acc_list.append(r_A_A_acc)
r_B_list.append(r_B)
r_cond.append(abs(r_A_A_acc[:]).max()) #x-x_acc Infinite norm of
r_B_cond.append(abs(r_B[:]).max())#b-Hx_acc Infinite norm of

print(' Dimension is 10,50,100 At the time of the x-x_acc The result of the calculation is ：\n', r_A_A_acc_list[9] ,r_A_A_acc_list[49] , r_A_A_acc_list[99])
print(' Dimension is 10,50,100 At the time of the b-Hx_acc The result of the calculation is ：\n',r_B_list[9], r_B_list[49], r_B_list[99])
print(' Dimension is 10,50,100 At the time of the x-x_acc The calculation result of infinite condition number of matrix ：\n', r_cond[9], r_cond[49], r_cond[99])
print(' Dimension is 10,50,100 At the time of the b-Hx_acc The calculation result of infinite condition number of matrix ：\n', r_B_cond[9], r_B_cond[49], r_B_cond[99])

``````

## 3. Output results

Enter what you want to calculate Hilbert The maximum dimension of a square matrix is OK , Leave the rest to the program .

`````` Please enter Hilbert The maximum dimension of a square matrix :100
Generating dimensions from 10,20……100 Of Hilbert The row norm condition number of a matrix ：
[35356847610517.12, 6.008376652086652e+18, 8.396589803249062e+18, 9.491653209312077e+19, 1.7763569870536153e+20, 1.9301974218850052e+21, 3.9847310708042826e+19, 1.3450693870678838e+20, 5.444272740462528e+19, 1.3244131088115743e+20]

``````

The output image here is as follows ：

Here is the output of the fourth question ：

`````` Dimension is 10,50,100 At the time of the x-x_acc The result of the calculation is ：
[[-2.54168641e-04]
[ 2.16242671e-03]
[-5.54656982e-03]
[ 5.08880615e-03]
[ 9.15527344e-04]
[-4.02832031e-03]
[ 1.46484375e-03]
[ 4.88281250e-04]
[-1.22070312e-04]
[-6.10351562e-05]] [[ 8.01768149e+02]
[ 3.33788188e+04]
[-1.59537467e+06]
[ 1.98594595e+07]
[-1.31128704e+08]
·················
[-4.04700000e+03]
[ 6.50000000e+01]
[-2.25000000e+01]
[ 3.30000000e+01]
[-7.30000000e+01]] [[ 1.05255071e+04]
[-1.31934071e+06]
[ 4.56146227e+07]
[-7.42843201e+08]
[ 6.95696228e+09]
[-4.13027099e+10]
················
[-4.79000000e+02]
[ 5.12100000e+03]
[-1.91900000e+03]
[ 1.77000000e+02]]
Dimension is 10,50,100 At the time of the b-Hx_acc The result of the calculation is ：
[[1.27400597e-05]
[2.15601768e-05]
[1.73587501e-05]
[1.43429989e-05]
[1.23056437e-05]
[1.08422262e-05]
[9.72981380e-06]
[8.84754702e-06]
[8.12566450e-06]
[7.52108032e-06]] [[ 13.49201973]
[  7.51301334]
[ -4.25849823]
[-14.49468816]
[-21.55733783]
[-26.46828937]
···············
[-28.3616173 ]
[-28.05340528]
[-27.75051982]
[-27.45291237]] [[-22.02594035]
[-22.62163018]
[-20.05848154]
[-17.70284588]
[-16.01064382]
[-14.79343569]
···············
[ -3.83066178]
[ -3.79759674]
[ -3.76500334]
[ -3.73286444]
[ -3.70119334]]
``````
`````` Dimension is 10,50,100 At the time of the x-x_acc The calculation result of infinite condition number of matrix ：
0.00554656982421875 7824513409.0 998313040247.0
Dimension is 10,50,100 At the time of the b-Hx_acc The calculation result of infinite condition number of matrix ：
2.1560176801216357e-05 38.404672581436365 22.62163018366762
``````

# Conclusion

Sharing is only for learning from each other , Discuss with each other , As the writing is wrong , Please forgive me a lot . Hope to learn from each other , Common progress , Welcome to leave a comment .