已知f(x)求x(f)的python代码

已知函数f(x),可以用插值法牛顿迭代法求x(f)。如果x和f(x)存在一对多的情况,x(f)就不是函数,要注意选取合适的域来避免这个情况。


插值法就是将输入的坐标用一个函数串起来。我们输入一系列的坐标(x,f(x)),那么x(f)的函数就是 scipy.interpolate.interp1d(x, f, kind="cubic")

举一个例子,已知\(a(t)=t^{2}+t\),求\(t(a)\)

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
def aoft(t):
    return t**2+t
t = np.linspace(0,2,10)
a =aoft(t)
tofa = interp1d(a, t, kind="cubic")
print(tofa(5),tofa(5)**2+tofa(5))

注意第六行,为了t(a)是函数,t取了大于0的情况。如果把第六行的取值范围改为-2到2,t(a)出来就是这样的:

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
def aoft(t):
    return t**2+t
t = np.linspace(-2,2,10)
a =aoft(t)
tofa = interp1d(a, t, kind="cubic")
aa=np.linspace(0,6,100)
plt.plot(a,t,"o",a,t,aa,tofa(aa),"--")
plt.xlabel(r'$a$')
plt.legend(['data','exact', 'interpolate'], loc='best')
plt.show()


牛顿迭代法解决的是已知f(x)=0,求x的问题。
首先我们需要猜一个解,\(x=x_{0}\)

然后,不停地计算\(x_{0}=x_0-\frac{f(x_0)}{f'(x_0)}\)
直到\(f(x_0)=0\)
这时\(x_{0}\)就是所求的解。
对于多根的方程,最开始猜的值不同,可能会得到不同的解。

例如,我需要解\(a(t)=t^{2}+t=5\)的解,对应的python代码是

def aoft(t):
    return t**2+t
def tofa(a):
    t = a**.5 #猜值
    def func(t):
        return aoft(t)-a
    for i in range(40):
        if(abs(func(t))<1e-6):break
        t = t-func(t)**2/(func(t+func(t))-func(t))
        if(i==39):print("for a: ",a,", the t_a finder does not converge\n")
    return t
print(tofa(5),tofa(5)**2+tofa(5))

如果想输入一系列的a得到对应的t,可以这样写

def tofan(a):
    def aoft(t):
        return t**2+t
    def tofa(a):
        t = -1 #猜值
        def func(t):
            return aoft(t)-a
        for i in range(40):
            if(abs(func(t))<1e-6):
                #print("for a: ",a,", the t_a : ",t," finder is convergent")
                break
            t = t-func(t)**2/(func(t+func(t))-func(t))
            if(i==39):
                print("for a: ",a,", the t_a: ",t,", finder does not converge!\n")
        return t
    if(type(a)==np.ndarray):
        tofa = [tofa(an) for an in a]
        return(np.array(tofa))
    else:
        return tofa(a)
a = np.linspace(0,10,1000)
plt.plot(a,tofan(a))

t的猜值大于0和小于0的图不一样,大家可以试试。
例子的函数很简单,我们直接画图比较一下

import matplotlib.pyplot as plt
import numpy as np
x=np.arange(-5,4,0.01)
y=np.arange(-1,10,0.01)# 转化为网格
x,y=np.meshgrid(x,y)
z=x**2+x-y
plt.contour(y,x,z,0)
Smilie Vote is loading.
✿ 阅读数:387  分类:其他

发表评论

电子邮件地址不会被公开。 必填项已用*标注

Captcha Code