給定區間記憶體在多根情況下的非線性方程求根

可愛的Y.y發表於2020-12-28
記憶體

Python 非線性方程求多根

非線性方程求根

常見的解決“給定區間的非線性方程求根”有二分法、牛頓法、割線法、不動點迭代法……
例如求 x 3 − x − 1 = 0 x^{3}-x-1=0 x3−x−1=0在區間 [ 1 , 2 ] [1,2] [1,2]內的根，我們就有多種方法進行求解。

二分法

#求x^3-x-1=0在區間[1,2]內的根
import numpy as np
from sympy import *
#二分法
def erfen(f,interval,tol):
    a=interval[0] #區間左側
    b=interval[1] #區間右側
    n=0
    while(b-a)/2>tol:
        c=(a+b)/2 #二分
        if f(c)==0:
            break
        if f(a)*f(c)<0: #替換b或者c的條件
            b=c
        else:
            a=c
        n+=1
    return c,n
f=lambda x:x**3-x-1
interval=[1,2]
tol=1e-5
print('二分法方程的解和迭代次數分別為：')
print(erfen(f,interval,tol))

牛頓法

import numpy as np
from sympy import *
#牛頓法
def newton(x0,f,df,tol):
    x=[]
    x_old=x0
    x.append(x0)
    error=1
    n=0
    while error>tol:
        x_new=x_old-f(x_old)/df(x_old)
        error=np.abs(x_new-x_old)
        x.append(x_new)
        x_old=x_new
        n+=1
    return x[-1],n
f=lambda x:x**3-x-1 
df=lambda x:3*x**2-1
print('牛頓法方程的解和迭代次數分別為：')
print(newton(1,f,df,tol))

割線法

import numpy as np
from sympy import *
#割線法
def gexian(x0,x1,f,tol):
    x=[]
    x.append(x0)
    x.append(x1)
    i=1
    error=1
    while error>tol:
        x_new=x[i]-f(x[i])*(x[i]-x[i-1])/(f(x[i])-f(x[i-1]))
        x.append(x_new)
        error=np.abs(x[i+1]-x[i])
        i+=1
    return x[-1],i
f=lambda x:x**3-x-1 
print('割線法方程的解和迭代次數分別為：')
print(gexian(1,2,f,tol))

但一般預設區間內僅有單根存在。
如果題目所給的方程在給定區間內未知根的個數，那又該如何進行計算呢？

非線性方程求多根

計算如下方程在區間[-10,10]中的所有根
∑ i = 1 10 k e − c o s k x s i n k x = 2 \sum_{i=1}^{10} ke^{-coskx}sinkx=2 i=1∑10​ke−coskxsinkx=2
這裡給出割線法求多根的一種方法，該方法適用於對於影像極為複雜難以看出零點的方程，我們可以運用該方法得到相應的零點個數，可以對複雜方程無需影像直接進行資料處理。

思路分析

首先考慮畫圖，從影像上直觀感受根的數量和區間之間的關係。
我們可以根據影像反映的各個根的區間對大區間進行相應的劃分變為小區間，在每個小區間內計算根的值，最後彙總得到大區間內的根的數量。
想進行更加準確的根區間的判斷，可以在上述根中求得兩根間的最小距離，並以此為區間劃分距離，再次對新劃分的小區間進行根的計算，判斷是否產生新根。若不產生新根，即說明已找到所有的根；若產生新根，則重複上述行動直至不產生新根。

程式碼

#上述問題程式
import numpy as np
import math
import matplotlib.pyplot as plt

#改寫方程為f(x)=Σk*np.exp(-np.cos(k*x))*np.sin(k*x)-2=0
def f(x):
    s=0
    for k in range(1,11):
        s=s+k*np.exp(-np.cos(k*x))*np.sin(k*x)
    return s-2

#畫影像初步看零點
plt.figure(figsize=(10, 10), dpi=70)
x=np.arange(-10.0,10.0,0.01)
plt.plot(x,f(x))
plt.axhline(y=0,ls="-",c="red") 
plt.show()


#精度
tol=1e-8
#零點個數
n=0
#儲存零點
y=[]

#以長度為1劃分[-10,10]區間
for j in range(20):
        #割線演算法gexian
    def gexian(x0,x1,f,tol):
        x=[]
        x.append(x0)
        x.append(x1)
        i=1
        error=1
        while error>tol:
            x_new=x[i]-f(x[i])*(x[i]-x[i-1])/(f(x[i])-f(x[i-1]))
            x.append(x_new)
            error=np.abs(x[i+1]-x[i])
            i+=1
            if x[i]>-10+1*(j+1) or x[i]<-10+1*j:
                break
            else:
                return x[-1]
    #判斷小區間內是否存在零點
    if gexian(-10+1*j,-10+1*(j+1),f,tol)!=None:
        y.append(gexian(-10+1*j,-10+1*(j+1),f,tol))
        n=n+1
print('零點的值為：',y)
print('零點的個數為：',n)

#求零點間的最小距離
k=len(y)
d=10
for i in range(k-1):
    dd=np.abs(y[i+1]-y[i])
    if dd<=d:
        d=dd
    else:
        d=d
print('零點間的最小距離為：',round(d,3))

#重新計算零點個數，以長度為d劃分[-10,10]區間
y=[]
n=0
geshu=int(20/d)
print('新的劃分割槽間個數為：',geshu)
for j in range(geshu+1):
        #割線演算法gexian
    def gexian(x0,x1,f,tol):
        x=[]
        x.append(x0)
        x.append(x1)
        i=1
        error=1
        while error>tol:
            x_new=x[i]-f(x[i])*(x[i]-x[i-1])/(f(x[i])-f(x[i-1]))
            x.append(x_new)
            error=np.abs(x[i+1]-x[i])
            i+=1
            if x[i]>-10+d*(j+1) or x[i]<-10+d*j:
                break
            else:
                return x[-1]
    #判斷小區間內是否存在零點
    if gexian(-10+d*j,-10+d*(j+1),f,tol)!=None:
        y.append(gexian(-10+d*j,-10+d*(j+1),f,tol))
        n=n+1
print('零點的值為：',y)
print('零點的個數為：',n)

#求零點間的最小距離
k=len(y)
d=10
for i in range(k-1):
    dd=np.abs(y[i+1]-y[i])
    if dd<=d:
        d=dd
    else:
        d=d
print('零點間的最小距離為：',round(d,3))

#重新計算零點個數，以長度為d劃分[-10,10]區間
y=[]
n=0
geshu=int(20/d)
print('新的劃分割槽間個數為：',geshu)
for j in range(geshu+1):
        #割線演算法gexian
    def gexian(x0,x1,f,tol):
        x=[]
        x.append(x0)
        x.append(x1)
        i=1
        error=1
        while error>tol:
            x_new=x[i]-f(x[i])*(x[i]-x[i-1])/(f(x[i])-f(x[i-1]))
            x.append(x_new)
            error=np.abs(x[i+1]-x[i])
            i+=1
            if x[i]>-10+d*(j+1) or x[i]<-10+d*j:
                break
            else:
                return x[-1]
    #判斷小區間內是否存在零點
    if gexian(-10+d*j,-10+d*(j+1),f,tol)!=None:
        y.append(gexian(-10+d*j,-10+d*(j+1),f,tol))
        n=n+1
print('零點的值為：',y)
print('零點的個數為：',n)

分析結果

上述程式碼的結果
我們可以看到在第二次迭代距離後根的個數就相應的穩定了下來，因此我們可以得到零點個數為79。

胡言亂語

PS：咳咳，我之後的小朋友們如果看到G老師的課還是用Python上的話，大家就會對這道題心知肚明知道我為什麼要用這種方法做了。造福小朋友吧嘻嘻~

相關文章