用Newton迭代,二分法,弦截法求根的通用程序
数值计算课程作业,供大家参考
逻辑思想参考:
https://wenku.baidu.com/view/cc3a5d0b0c22590103029d06.html
语言:python
程序介绍:本程序适用于一元方程数值求解,要求一元方程对应的函数在定义域内连续,但不要求可导。本程序可以处理方程有一个根,两个根以及多个根的情形,程序自动寻找x0的初值,但如果函数在x轴上下波动特别剧烈,程序可能无法找到全部根,这时需要手动增加分辨率。
import math
import numpy as np
import matplotlib.pyplot as plt
import sympy
def Findx1(a,b): #a b 是定义域,这个函数用于寻找x1作为迭代的初始值
h = 100.
global dx
dx = (b-a)/h
x = np.arange(a,b+dx,dx)
f1 = (abs(Fun(x+dx))-abs(Fun(x)))/dx #取绝对值,让函数图形位于x轴上方
f2 = (abs(Fun(x))-abs(Fun(x-dx)))/dx
ff = f1*f2 #负数代表函数在x位置有极值
mask = ff<0
x0 = x[mask]
mask2 = (abs(Fun(x0-dx))-abs(Fun(x0)))>0 #如果函数绝对值递增出现极值,肯定不是解,舍去
x1 = x0[mask2]
print(x1)
return x1
def Newton(xk): #牛顿法计算精确值,每次只计算一个根
dfx1 = abs(Fun(xk+dx)-Fun(xk))/dx #前面几行比较x0 和x0+h哪个更容易收敛,取导数大的那个
dfx2 = abs(Fun(xk+2.*dx)-Fun(xk+dx))/dx
if dfx1<dfx2:
xk = xk+dx
xk1 = xk+dx
i = 0 #
dx0 = 1. # dx0用于求导,这里是初始值。如果已知函数导数表达式,这里可以替换掉
while(abs(dx0)>10e-7 and i < 10): #控制循环次数
dfx0 = Dfun(xk) #x点处的导数
xk1 = xk-Fun(xk)/dfx0 #计算新的x
dx0 = xk1-xk #新的dx,用于求导
xk = xk1
i = i+1
if abs(Fun(xk))>10e-2: #判断f(x)是不是在x轴上,如果不是返回99999
return 999999
return xk
def Secant(xk): #弦截法,基本上和牛顿法一样
dfx1 = abs(Fun(xk+dx)-Fun(xk))/dx #
dfx2 = abs(Fun(xk+2.*dx)-Fun(xk+dx))/dx
if dfx1<dfx2:
xk = xk+dx
xk1 = xk+dx
i = 0 #
dx0 = xk1-xk # dx0用于求导,这里是初始值
while(abs(dx0)>10e-7 and i < 10): #控制循环次数
dfx0 = (Fun(xk+dx0)-Fun(xk))/dx0 #x点处的导数
xk1 = xk-Fun(xk)/dfx0 #计算新的x
dx0 = xk1-xk #新的dx,用于求导 ,这行和牛顿法略有不同
xk = xk1
i = i+1
if abs(Fun(xk))>10e-4: #判断f(x)是不是根,如果不是返回99999,然后主程序里面可以将这个值过滤掉
return 999999
return xk
def bisection(xk): #二分法。如果函数二阶导数为零,则牛顿法失效,换二分法
xk1 = xk+dx
err = 1.
k = 1
while (err>10e-15 and k<20):
if abs(Fun(xk))>abs(Fun(xk1)):
xk = (xk+xk1)/2.
else:
xk1 = (xk+xk1)/2.
k = k+1
err = xk1-xk
return xk
def Printresult(a): # 这部分先判断用哪种方法解方程, 然后调用解方程的函数,过滤无效解
ddfx1 = (Fun(a+dx)-2.*Fun(a)+Fun(a-dx))/dx/dx #a点处的二阶差商
ddfx2 = (Fun(a+2.*dx)-2.*Fun(a+dx)+Fun(a))/dx/dx #a+dx点处的二阶差商
mask = ddfx1*ddfx2<0 #判断二阶导数的变化,如果变号用二分法,否则用牛顿法
b = a[mask] #在这个根处二阶导数变号,用二分法求根
mask2 = ddfx1*ddfx2>= 0 #不变号的用牛顿法
a2 = a[mask2]
print(mask2)
for i in range(0,a2.size):
xx = Newton(a2[i])
if abs(xx-999999)<1:
continue
print('Newton',xx)
for j in range(0,b.size):
yy = bisection(b[j])
if abs(yy-999999)<1:
continue
print('bisection',yy)
return
def Dfun(x):
df=x*x-1.
return df
def Fun(x): #函数形式
f = x**3/3.-x#+50.*np.sin(x)+24.5 #
return f
a = Findx1(-10.,10.) #括号内是定义域,找到根附近的值,作为初值x0
Printresult(a)
xxx = np.arange(-10,10,0.1) #最后画函数图,检验解是否和函数图一致
plt.figure(1)
plt.plot(xxx,Fun(xxx))
plt.show()