突然想起来可以做做数值分析的作业，于是把室友的数值分析作业拿过来练手，写一篇博客分享一下。其实我这个菜鸟把程序写复杂了，有很多可以简化的地方，由于本菜鸟其它作业还没写完，就不去简化了，大家可以自行改正啦。

---

数值分析上机题

首先说一下自己的疑惑，对于第一题python怎么实现对ln(x)直接调用求导呢？是直接用泰勒展开后对多项式求导吗？第二题是用Newton迭代法，要求出迭代初始值在什么范围可以得到收敛解，这里如果用迭代程序去实现的话也会比较麻烦，那有没有更好的方法去求解呢？希望有大神可以帮忙留言解惑，多谢。

一、题1——对数的近似求解

1.题目描述

**题目：**这里偷下懒，直接贴出原题的截图吧。

2.python实现

不多BB，程序在这：

import numpy as np
from sympy import * #用于求导积分等科学计算
import math as mx = Symbol('x')#x 变量
t = Symbol('t')
y1 = 1/(1+x) #公式
y2 = -1/(1+x) #公式
y3 = 2/(1-x**2) #公式def func(m):res = mfor j in range(1,m):res *= jreturn resdef 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的值：', m.log(2, m.e))
print('方程ln(1+x)的10阶泰勒展开计算结果为：', A,'\n','估计误差为：', abs(m.log(2, m.e)-A))
print('方程ln(1/(1+x))的10阶泰勒展开计算结果为：', B,'\n','估计误差为：', abs(m.log(2, m.e)-B))
print('方程ln((1+x)/(1-x))的10阶泰勒展开计算结果为：', C,'\n','估计误差为：', abs(m.log(2, m.e)-C))

3.输出结果

ln2的值： 0.6931471805599453
方程ln(1+x)的10阶泰勒展开计算结果为： 0.645634920634921 估计误差为： 0.0475122599250246
方程ln(1/(1+x))的10阶泰勒展开计算结果为： 0.693064856150793 估计误差为： 8.23244091517905e-5
方程ln((1+x)/(1-x))的10阶泰勒展开计算结果为： 0.693146047390827 估计误差为： 1.13316911820593e-6

一、题2——方程求根的Newton法

1.题目描述

**题目：**还是截图方便。

2.python实现

这个程序写的不好，由于写完后，exp（）函数老是报错说：整型数据不是exp对象的属性，改了之后也没实现自己的想法，那就这样吧，没时间啦，大家自己搞吧。。。。。

import numpy as np
from sympy import * #用于求导积分等科学计算
import math as mdef f(x):return 3*x**2 - m.exp(x) #该方程有3个根，初步估计在[-1,0],[0,1],[3,4]def fdiff(x):return 6*x - m.exp(x)def g(x):return 6*x - m.exp(x)#该方程有3个根，初步估计在[0,1],[2,3]def gdiff(x):return 6 - m.exp(x)a = float(input('请输入计算区间的下界a(浮点型）: '))
b = float(input('请输入计算区间的上界b(浮点型）: '))
c = float(input('请输入迭代初始值(浮点型）: '))
n = input('请输入要求解的函数，f代表f(x),g代表g(x): ')if n =='f':if f(a) * f(b)> 0:print('在此区间内函数有多个根或者无根')elif f(a) * f(b) == 0:print('f(a) = ', '%f'%f(a))print('f(b) = ', '%f'%f(b))else:fcount = 0y = c - f(c) / fdiff(c)while (abs(c - y) >= 0.5e-9) & (fdiff(c) != 0):x2 = c - f(c) / fdiff(c)y = cc = x2fcount += 1print('函数给出的求根区间是：', [a, b])print("函数f(x）的Newton迭代次数：%f,函数f(x）的迭代计算所得的根为:%.8f"%(fcount,c))elif n =='g':if g(a) * g(b)> 0:print('在此区间内函数有多个根或者无根')elif g(a) * g(b) == 0:print('g(a) = ', '%f'%g(a))print('g(b) = ', '%f'%g(b))else:gcount = 0y = c - g(c) / gdiff(c)while (abs(c - y) >= 0.5e-9) & (gdiff(c) != 0):x2 = c - g(c) / gdiff(c)y = cc = x2gcount += 1print('函数给出的求根区间是：', [a, b])print("函数g(x）的Newton迭代次数：%f,函数g(x）的迭代计算所得的根为:%.8f"%(gcount,c))

3.输出结果

这里说明一下，输入的[a, b]是你想判断该区间有没有根；c是迭代初始值；f代表f（x)，g代表g(x)；这几个参数请自己输入。

请输入计算区间的下界a(浮点型）: -1.0
请输入计算区间的上界b(浮点型）: 4.0
请输入迭代初始值(浮点型）: -3.0
请输入要求解的函数，f代表f(x),g代表g(x): f
函数给出的求根区间是： [-1.0, 4.0]
函数f(x）的Newton迭代次数：7.000000,函数f(x）的迭代计算所得的根为:-0.45896227

一、题3——Hilbert矩阵的条件数

1.题目描述

**题目：**你懂的。

2.python实现

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['simhei']n = int(input('请输入Hilbert方阵的最大维数:' ))
def max_sum_rows(X):sum_row_list1 = []for i in range(len(X)):count = 0for 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 * a2Hilbert_cond.append(H_cond)#计算10,20……100的无穷范数条件数
Hilbert_cond_test = []
for j in range(n):if (j+1)%10 == 0:Hilbert_cond_test.append(Hilbert_cond[j])
print('生成维数从10,20……100的Hilbert矩阵的行范数条件数：\n', Hilbert_cond_test)plt.title('Hilbert矩阵维数与条件数对数之间的关系')
plt.plot((list(range(1,n+1))), np.log(Hilbert_cond),c ='b', marker='*',label='拟合曲线')
plt.legend()
plt.xlabel('矩阵维度n')
plt.xticks(np.arange(0, n+1, 5))
plt.yticks(np.arange(0, 55, 5))
plt.ylabel('log(cond)')
plt.show()#求解Hilbert方程的解和对应的无穷条件数
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))*1X = 1. / (np.arange(1, i + 1) + np.arange(0, i)[:, np.newaxis])B = X@AA_acc = np.linalg.inv(X)@Br_A_A_acc = A - A_accr_B = B - X @ A_accr_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的无穷范数r_B_cond.append(abs(r_B[:]).max())#b-Hx_acc的无穷范数print('维数为10,50,100时的x-x_acc计算结果：\n', r_A_A_acc_list[9] ,r_A_A_acc_list[49] , r_A_A_acc_list[99])
print('维数为10,50,100时的b-Hx_acc计算结果：\n',r_B_list[9], r_B_list[49], r_B_list[99])
print('维数为10,50,100时的x-x_acc矩阵的无穷条件数计算结果：\n', r_cond[9], r_cond[49], r_cond[99])
print('维数为10,50,100时的b-Hx_acc矩阵的无穷条件数计算结果：\n', r_B_cond[9], r_B_cond[49], r_B_cond[99])

3.输出结果

输入你想计算的Hilbert方阵的最大维数就行，其它交给程序。

请输入Hilbert方阵的最大维数:100
生成维数从10,20……100的Hilbert矩阵的行范数条件数：[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]

这里输出的图片如下：

这里是第四问的输出结果：

维数为10,50,100时的x-x_acc计算结果：[[-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]]维数为10,50,100时的b-Hx_acc计算结果：[[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]]

维数为10,50,100时的x-x_acc矩阵的无穷条件数计算结果：0.00554656982421875 7824513409.0 998313040247.0
维数为10,50,100时的b-Hx_acc矩阵的无穷条件数计算结果：2.1560176801216357e-05 38.404672581436365 22.62163018366762

结语

分享出来仅供相互学习，相互探讨，如所写有误，请多多包涵。希望能相互学习，共同进步，欢迎各位大佬留言评论。

shiyou的数值分析作业相关推荐

1. 欧拉编程c语言作业数值分析,数值分析作业 欧拉 龙格库塔

   <数值分析作业 欧拉 龙格库塔>由会员分享,可在线阅读,更多相关<数值分析作业 欧拉 龙格库塔(5页珍藏版)>请在人人文库网上搜索. 1.9.2对初值问题试用欧拉方法取步长h= ...

2. 龙贝格数值分析作业c语言,数值分析龙贝格实验报告.doc

   数值分析龙贝格实验报告 实验三 龙贝格方法 [实验类型] 验证性 [实验学时] 2学时 [实验内容] 1.理解龙贝格方法的基本思路 2.用龙贝格方法设计算法,编程求解一个数值积分的问题. [实验前的预 ...

3. 求f(x)=根号x的三次样条插值和分段线性插值c语言代码,【数值分析|三次样条插值法《数值分析》上机实验作业】...

   『易坊知识库摘要_数值分析|三次样条插值法<数值分析>上机实验作业』2.要求(1)满足自然边界条件s?(0.2)?s?(1.0)?0:(2)满足第一类边界条件s?(0.2)?0.20271 ...

4. 艺术外包兼职网站app推荐_做兼职的艺术

   艺术外包兼职网站app推荐 Hi everyone, I'm going to share some of my past experiences doing part-time work, as w ...

5. matlab迭代求解,[基于matlab平台的三种迭代法求解矩阵方程]matlab迭代法求方程的根...

    数值分析第二次作业 学院:电子工程学院 基于matlab平台的三种迭代法求解矩阵方程组 求解系数矩阵由16阶Hilbert方程组构成的线性方程组的解,其中右端项为[2877/851,3491/14 ...

6. matlab 高斯迭代法求解,高斯迭代法matlab算例

   Matlab 线性方程组的迭代解法 Gauss-Seidel 迭代法 Matlab 线性方程组的迭代解法 Gauss-Seidel 迭代法实验报告 1.熟悉 Gauss-Seidel 迭代法,并编写 ...

7. 研究生日迹-201710月

   研究生日迹-201710月 20171004 中秋 这个国庆节在学校玩的很嗨,吃了很多好吃的,还去万达电玩城玩,看了电影,游戏是打足了,歌也唱过了,总之很尽兴很开心. 不过这几天钱花的过多,已经用了6 ...

8. 数值分析matlab实验报告,数值分析第一次作业matlab实验报告.doc

   数值分析第一次作业matlab实验报告.doc 几种线性方程组迭代算法的MATLAB实现和性能比较用有限差分方法(五点差分格式)求解正方形域上的Poisson方程边值问题用MATLAB语言编写算法程序 ...

9. 数值分析实习作业（各种插值函数与积分公式的python代码实现）

   数值分析实习作业 1.对函数 f(x) f ( x ) f(x)进行插值 f(x)=11+x2 f ( x ) = 1 1 + x 2 f(x)=\frac{1}{1+x^{2}} (1)令插值节点为 ...

10. 数值分析的matlab答案,Matlab作业3(数值分析)答案.doc

    Matlab作业3(数值分析)答案 Matlab作业3(数值分析) 机电工程学院 (院.系) 专业 班 组 学号 姓名 实验日期 教师评定 计算多项式乘法(x2+2x+2)(x2+5x+4). 答: ...

最新文章

1. 介绍 Saltstack批量管理文件和计划任务
2. Xamarin.Android 使用 Encoding.GetEncoding(GB2312) 报错解决方案
3. 【深度学习】Coursera的TensorFlow课程练习题精华部分
4. python批量做线性规划（每次的约束条件参数有变换）
5. 实用脚本！Python 提取 PDF 指定内容生成新文件！
6. oracle内外链接混合用,混合在一起通过连接，内部连接和总结与Oracle
7. centos7 安装webmin
8. C语言 · 数的读法
9. Dubbo(七)使用SpringBoot搭建dubbo消费者工程
10. c语言语句的使用形式,C语言如何使用print语句
11. Oracle网络配置用到的sqlnet.ora,tnsnames.ora,listener.ora文件
12. 异型烟分拣 机器人_细支烟、标准烟共线分拣可行性研究
13. 使用脚本解决fstab挂载失败不能正常启动问题
14. java代码生成springdao_可一键生成dao、表、controller等几十种的代码生成器源码分享...
15. 数据安全--3--数据安全5A之授权
16. 戴建业老师对李白和杜甫的讨论
17. xbox360链接pc_如何在Windows PC上使用Xbox 360控制器
18. 什么是interop
19. 视频教程-Spring Boot实战入门视频课程-Java
20. 少年歌行游戏一直显示连接服务器,少年歌行第二季无法上线的原因找到了，不是没做好，而是限制太多...

热门文章

1. SPSS学习（五）独立样本t检验
2. 【gflags】【gflags实践】【gflags的学习使用记录】
3. 三轴加速度传感器和六轴惯性传感器_六轴传感器和三轴传感器的区别
4. STM32F407 USB CDC调试与经验总结
5. 极路由第三方插件大全_极路由极硬货HC5663春节折腾记
6. ①编写一个程序，实现文件的复制。②写一个加密程序，对文件1.txt进行加密。它从输入流中读入一个密钥，并将文件密文输出。③ 编写一个应用程序,向数据文件out.dat中输入100个1000以内的随机整
7. svn执行reflash/cleanup报错wc.db解决办法
8. matlab常用插值函数
9. 华为鸿蒙deveco studio编译时提示Browserslist: caniuse-lite is outdated的解决办法
10. 四川途志:短视频营销公司做视频广告投放有技巧吗？