Мой код не дает мне правильное значение - PullRequest
0 голосов
/ 14 февраля 2019

Я делаю это для домашней работы (извините за мой английский). Проблема в жирных строках: точно проблема с idk, но похоже, что xi + 1 принимает неправильные значения (я объясню в конце):

import numpy as np
from astropy.table import Table
def f(x): #function
    return x**3+2*x**2+10*x-20
#Initial Values, I put matrix of 100 because this method is for aproximation so we don't know how many iterations it will take
i=0
l=0.00004
e=np.zeros([100,1])
x=np.zeros([100,1])
a=np.zeros([100,1])
b=np.zeros([100,1])
G=np.zeros([100,1])
S=np.zeros([100,1])
a[0]=1
b[0]=3
    #Aproximation method
def FP(x,i,a,b):
    G[i]=f(b[i])
    F[i]=f(a[i])
    x[i]=(a[i]*G[i]-b[i]*F[i])/(G[i]-F[i])
    e[i]=abs(x[i]-a[i])
    while e[i]>l:
        if f(a[i])*f(x[i])>0:
            a[i+1]=x[i]
            b[i+1]=b[i]
            G[i+1]=f(b[i+1])
            F[i+1]=f(a[i+1])
#HERE
            x[i+1]=(a[i+1]*G[i+1]-b[i+1]*F[i+1])/(G[i+1]-F[i+1])
            if f(x[i])*f(x[i+1])>0:
                G[i+1]=G[i+1]/2
            e[i+1]=abs(x[i+1]-x[i])
            i=i+1
        elif f(a[i])*f(x[i])<0:
            a[i+1]=a[i]
            b[i+1]=x[i]
            G[i+1]=f(b[i+1])
            F[i+1]=f(a[i+1])
#AND HERE
            x[i+1]=(a[i+1]*G[i+1]-b[i+1]*F[i+1])/(G[i+1]-F[i+1])
            if f(x[i])*f(x[i+1])>0:
                F[i+1]=F[i+1]/2
            e[i+1]=abs(x[i+1]-x[i])
            i=i+1
        else:
            break
 #This is for the table so ignore it     
    im=range(0,i+1)
    xm=np.around(x[0:i+1],4)
    am=np.around(a[0:i+1],4)
    bm=np.around(b[0:i+1],4)
    em=np.around(e[0:i+1],4)
    Gm=np.around(G[0:i+1],4)
    Fm=np.around(F[0:i+1],4)
    fxm=np.around(f(x[0:i+1]),4)
    t = Table([im,am,bm,Fm,Gm,xm,fxm,em], names=('im','ai', 'bi', 'F','G','xi+1','f(xi+1)','e'))
    print(t) 

И когда я оцениваю FP (x, i, a, b) во второй итерации, что-то не так

im ai [1] bi [1]  F [1]  G [1] xi+1 [1] f(xi+1) [1] e [1] 
--- ------ ------ ------- ----- -------- ----------- ------
  0    1.0    3.0    -7.0  55.0   1.2258     -2.8948 0.2258
  1 1.2258    3.0 -2.8948  27.5   1.3145     -1.1275 0.0887

Xi + 1должно быть (ai Gi-bi Fi) / (Gi + Fi) (должно принимать значения одной строки), поэтому xi + 1 = 1,39478 (приблизительно), но программа дает мне 1,3145 (приблизительно),Я знаю ответ, потому что я сделал то же самое в техасском n-spire CX CAS (там действительно легче), и мой учитель получил тот же ответ.Спасибо

Это метод, который я использую Алгоритм рис Fi = f (ai) и Gi = f (bi) Theполный код на испанском (если вам нужно запустить)

import numpy as np
from astropy.table import Table
def f(x): #función 
    return x**3+2*x**2+10*x-20
#Valores iniciales
i=0
l=0.00004
e=np.zeros([100,1])
x=np.zeros([100,1])
a=np.zeros([100,1])
b=np.zeros([100,1])
G=np.zeros([100,1])
S=np.zeros([100,1])
a[0]=1
b[0]=3

#Aproximación por método de posición falsa
def FP(x,i,a,b):
    G[i]=f(b[i])
    F[i]=f(a[i])
    x[i]=(a[i]*G[i]-b[i]*F[i])/(G[i]-F[i])
    e[i]=abs(x[i]-a[i])
    while e[i]>l:
        if f(a[i])*f(x[i])>0:
            a[i+1]=x[i]
            b[i+1]=b[i]
            G[i+1]=f(b[i+1])
            F[i+1]=f(a[i+1])
            x[i+1]=(a[i+1]*G[i+1]-b[i+1]*F[i+1])/(G[i+1]-F[i+1])
            if f(x[i])*f(x[i+1])>0:
                G[i+1]=G[i+1]/2
            e[i+1]=abs(x[i+1]-x[i])
            i=i+1
        elif f(a[i])*f(x[i])<0:
            a[i+1]=a[i]
            b[i+1]=x[i]
            G[i+1]=f(b[i+1])
            F[i+1]=f(a[i+1])
            x[i+1]=(a[i+1]*G[i+1]-b[i+1]*F[i+1])/(G[i+1]-F[i+1])
            if f(x[i])*f(x[i+1])>0:
                F[i+1]=F[i+1]/2
            e[i+1]=abs(x[i+1]-x[i])
            i=i+1
        else:
            break
    #Esto es para la tabla        
    im=range(0,i+1)
    xm=np.around(x[0:i+1],4)
    am=np.around(a[0:i+1],4)
    bm=np.around(b[0:i+1],4)
    em=np.around(e[0:i+1],4)
    Gm=np.around(G[0:i+1],4)
    Fm=np.around(F[0:i+1],4)
    fxm=np.around(f(x[0:i+1]),4)
    t = Table([im,am,bm,Fm,Gm,xm,fxm,em], names=('i','ai', 'bi', 'F','G','xi+1','f(xi+1)','e'))
    print(t) 
FP(x,i,a,b)
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...