一、目的与要求
(一)目的通过设计、编制、调试2~3个用数值方法求方程根的程序,加深对方程求根的数值计算方法及有关的基础理论知识的理解。
(二)要求 用编程语言实现二分法、Newton迭代法、弦截法求方程根的程序。
二、示例
1、问题已知f(a)与f(b)异号,用二分法求方程f(x)=0所在的根。
2、程序中变量说明 (略)
3、源程序清单及运行结果(略)
4、按以上4点要求编写上机实验报告。
做法:
1、用二分法求的根。
2、用牛顿(Newton)迭代法求在附近的根。
3、用弦截法求的根。
代码:
import numpy as np
import matplotlib.pyplot as plt
# from scipy.interpolate import spline
N = 10
TOL = 1e-5
# 模拟函数为func(x)
def func(x):
return x**3 - x - 1
# 导函数为f(x)
def f(x):
return 3*(x**2) - 1
#二分法
def Bindary():
low,high = 1,2
for i in range(N):
g = (low + high) / 2
if func(g) == 0 or (high-low)/2 < TOL:
return g
elif func(low)*func(g) > 0:
low = g
else:
high = g
# 返回二分法第N次计算所得的根
return g
#牛顿法
def Newton():
low, high = 1, 2
NewtonApprox = np.zeros(N+1)
NewtonApprox[0] = low
if func(low)==0:
return low
for i in range(1,N+1):
NewtonApprox[i] = NewtonApprox[i-1] - func(NewtonApprox[i-1])/f(NewtonApprox[i-1])
# 返回牛顿法N次计算所得的根
return NewtonApprox
#弦截法
def Secant():
low, high = 1, 2
SecantApprox = np.zeros(N+1)
SecantApprox[0] = low
SecantApprox[1] = low+0.1
if func(low)==0:
return low
for i in range(2,N+1):
SecantApprox[i] = SecantApprox[i-1] - func(SecantApprox[i-1])*(SecantApprox[i-1]-SecantApprox[i-2])/(func(SecantApprox[i-1])-func(SecantApprox[i-2]))
return SecantApprox
#画图观察各次求出的根的分布情况
def Draw():
Bind = Bindary()
print("二分法求解方程的近似根为 ",Bind)
#得到牛顿法 弦截法结果数组
Newton1 = Newton()
Secant1 = Secant()
x = 2*np.random.rand(10)
#思路:先在图像上画出原函数直线 y=0水平直线 再将牛顿法 弦截法结果数组各点画到图像上
x1 = np.linspace(1,2,502)
y1 = x1**3 - x1 - 1
plt.plot(x1,y1)
plt.xlabel('X',size=15,horizontalalignment='right')
plt.ylabel('Y',size=15,verticalalignment='top')
plt.axhline(y=0,color='b',linestyle='--')
plt.annotate(r'y = 0', xy=(2,0), xytext=(1.7, 0.4),size=15,
arrowprops=dict(facecolor='black', shrink=0.2, width=0.6))
plt.annotate(r'y = x**3 - x - 1', xy=(1.8,3), xytext=(1.3, 3.4),size=15,
arrowprops=dict(facecolor='black', shrink=0.2, width=0.6))
plt.scatter(Newton1[0:N],func(Newton1[0:N]),color='black')
plt.scatter(Secant1[0:N-1], func(Secant1[0:N-1]), color='orange')
plt.plot(1.1, 4.0, color="black", linewidth=2.5, linestyle='solid', label="Newton")
plt.plot(1.1, 3.5, color="orange", linewidth=2.5, linestyle='solid', label="Secant")
plt.legend(loc='upper left')
plt.title('Antiderivative graph')
plt.savefig('Extract a root',dpi=600)
plt.show()
if __name__ == '__main__':
Draw()
程序缺点:
复用性有些差