当前位置:网站首页>Shiyou's numerical analysis assignment
Shiyou's numerical analysis assignment
2020-11-08 08:04:00 【osc_4x0ulctb】
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 .
List of articles
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 .
版权声明
本文为[osc_4x0ulctb]所创,转载请带上原文链接,感谢
边栏推荐
- 1.深入Istio:Sidecar自动注入如何实现的?
- Littlest JupyterHub| 02 使用nbgitpuller分发共享文件
- M 端软件产品设计思虑札记 - 知乎
- Qt混合Python开发技术:Python介绍、混合过程和Demo
- Android Basics - RadioButton (radio button)
- Speed up your website with jsdelivr
- On the stock trading of leetcode
- Abnormal + Abstract
- iOS上传App Store报错:this action cannot be completed -22421 解决方案
- Cryptography - Shangsi Valley
猜你喜欢
C++基础知识篇:C++ 基本语法
Do you really understand the high concurrency?
wanxin金融
PCR and PTS calculation and inverse operation in TS stream
SQL Server 2008R2 18456错误解决方案
More than 50 object detection datasets from different industries
Face recognition: attack types and anti spoofing techniques
来自不同行业领域的50多个对象检测数据集
麦格理银行借助DataStax Enterprise (DSE) 驱动数字化转型
Cloud Alibabab笔记问世,全网详解仅此一份手慢无
随机推荐
技术人员该如何接手一个复杂的系统?
GET,POST,PUT,DELETE,OPTIONS用法与说明
0.计算机简史
Goland 编写含有template的程序
C/C++编程笔记:C语言相比其他编程语言,有什么不一样的优势?
C语言I博客作业03
1. In depth istio: how is sidecar auto injection realized?
C language I blog assignment 03
Six key points of data science interview
搜索引擎的日常挑战_4_外部异构资源 - 知乎
nvm
Lay UI left tree Dtree right list table
Judging whether paths intersect or not by leetcode
归纳一些比较好用的函数
5g + Ar out of the circle, China Mobile Migu becomes the whole process strategic partner of the 33rd China Film Golden Rooster Award
Golang anonymous structure member, named structure member, inheritance, composition
GoLand writes a program with template
洞察——风格注意力网络(SANet)在任意风格迁移中的应用
c# 表达式树(一)
高并发,你真的理解透彻了吗?