在Python中打印Bulirsch-Stoer算法的解决方案

问题描述 投票：0回答：1

在下面的代码中，为什么我不断出现此错误：

TypeError：“ NoneType”对象不可下标

对于我的PrintSoln(X, Y, freq)函数？我该如何解决？

我正在尝试解决IVP：2*(y'') + y' + (y/0.45) = 9，初始条件为y = 0 and y' = 0。'''

from numpy import zeros, float, sum
import math

def integrate_BulirschStoeir(F, x, y, xStop, tol):
    def midpoint(F, x, y, xStop, nSteps):
        h = (xStop - x)/ nSteps
        y0 = y
        y1 = y0 + h*F(x, y0)
        for i in range(nSteps-1):
            x = x + h
            y2 = y0 + 2.0*h*F(x, y1)

            y0 = y1
            y1 = y2

        return 0.5*(y1 + y0 + h*F(x, y2))

    def richardson(r, k):
        for j in range(k-1,0,-1):
            const = (k/(k - 1.0))**(2.0*(k-j))
            r[j] = (const*r[j+1] - r[j])/(const - 1.0)
        return     
        kMax = 51
        n = len(y)
        r = zeros((kMax, n), dtype=float)

        # Start with two integration steps
        nSteps = 2
        r[1] = midpoint(F, x, y, xStop, nSteps)
        r_old = r[1].copy()

        # Increase the number of integration points by 2 and refine result by Richardson extrapolation
        for k in range(2,kMax):
            nSteps = 2*k
            r[k] = midpoint(F, x, y, xStop, nSteps)
            richardson(r,k)

            # Compute RMS change in solution
            e = sqrt(sum((r[1] - r_old)**2)/n)

            # Check for convergence
            if e < tol: 
                return r[1]
            r_old = r[1].copy()
        print("Midpoint method did not converge")

# Bulirsch-Stoer Algorithm:- 

''' X, Y = bulStoer(F, x, y, xStop, H, tol=1.0e-6).
    Simplified Bulirsch-Stoer method for solving the
    initial value problem {y}’ = {F(x,{y})}, where {y} = {y[0],y[1],...y[n-1]}
    x, y = initial conditions
    xStop = terminal value of x
    H = increment of x at which results are stored
    F = user-supplied function that returns the array F(x,y) = {y’[0],y’[1],...,y’[n-1]} '''

from numpy import array

def bulStoer(F, x, y, xStop, H, tol=1.0e-6):
    X = []
    Y = []
    X.append(x)
    Y.append(y)
    while x < xStop:
        H = min(H,xStop - x)
        y = integrate_BulirschStoeir(F, x, y, x + H, tol)   # Midpoint method
        x = x + H
        X.append(x)
        Y.append(y)
    return array(X), array(Y)

def printSoln(X, Y, freq):
        def printHead(n):
            print("\n x ")
            for i in range (n):
                print("y[",i,"]")
            print
        def printLine(x, y, n):
            print("%13.4e"% x)
            for i in range (n):
                print("%13.4e"% y[i])
            print
        m = len(Y)
        try: n = len(Y[0])
        except TypeError: n = 1
        if freq == 0: freq = m
        printHead(n)
        for i in range(0,m,freq):
            printLine(X[i], Y[i], n)
        if i != m - 1: printLine(X[m - 1], Y[m - 1], n)

# Code:-

from numpy import array, zeros

def F(x, y):
    F = zeros(2)
    F[0] = y[1]
    F[1] =(-y[1] - (y[0]/0.45) + 9.0)/2.0
    return F

x = 0.0
xStop = 10.0
H = 0.5
y = array([0.0, 0.0])
freq = 1
X, Y = bulStoer(F, x, y, xStop, H)
printSoln(X, Y, freq)
python algorithm debugging ode
1个回答
-1
投票

  1. “ def integration_BulirschStoeir（F，x，y，xStop，tol）”（从“ kMax = 51”到函数结尾）存在缩进问题。因此，在“ X，Y = bulStoer（F，x，y，xStop，H）”和Y bacame的矩阵之后，Y中没有新数据，该矩阵的第一行为[0,0]，其他行为None。

  2. [[此外，“ def printLine”中的打印在打印功能中未添加“ end =''”导致了列向量，并添加了“ end =''”以使结果更具可读性。

    from numpy import zeros, float, sum, sqrt import math def integrate_BulirschStoeir(F, x, y, xStop, tol): def midpoint(F, x, y, xStop, nSteps): h = (xStop - x)/ nSteps y0 = y y1 = y0 + h*F(x, y0) for i in range(nSteps-1): x = x + h y2 = y0 + 2.0*h*F(x, y1) y0 = y1 y1 = y2 return 0.5*(y1 + y0 + h*F(x, y2)) def richardson(r, k): for j in range(k-1,0,-1): const = (k/(k - 1.0))**(2.0*(k-j)) r[j] = (const*r[j+1] - r[j])/(const - 1.0) return kMax = 51 ## Indentation problem from here n = len(y) r = zeros((kMax, n), dtype=float) # Start with two integration steps nSteps = 2 r[1] = midpoint(F, x, y, xStop, nSteps) r_old = r[1].copy() # Increase the number of integration points by 2 and refine result by Richardson extrapolation for k in range(2,kMax): nSteps = 2*k r[k] = midpoint(F, x, y, xStop, nSteps) richardson(r,k) # Compute RMS change in solution e = sqrt(sum((r[1] - r_old)**2)/n) # Check for convergence if e < tol: return r[1] r_old = r[1].copy() print("Midpoint method did not converge") # Bulirsch-Stoer Algorithm:- ''' X, Y = bulStoer(F, x, y, xStop, H, tol=1.0e-6). Simplified Bulirsch-Stoer method for solving the initial value problem {y}’ = {F(x,{y})}, where {y} = {y[0],y[1],...y[n-1]} x, y = initial conditions xStop = terminal value of x H = increment of x at which results are stored F = user-supplied function that returns the array F(x,y) = {y’[0],y’[1],...,y’[n-1]} ''' from numpy import array def bulStoer(F, x, y, xStop, H, tol=1.0e-6): X = [] Y = [] X.append(x) Y.append(y) while x < xStop: H = min(H,xStop - x) y = integrate_BulirschStoeir(F, x, y, x + H, tol) # Midpoint method x = x + H X.append(x) Y.append(y) return array(X), array(Y) def printSoln(X, Y, freq): def printHead(n): print("\n x ",end=" ") ## end=" " here for i in range(n): print(" y[",i,"] ",end=" ") ## end=" " here print() def printLine(x, y, n): print('{:13.4e}'.format(x), end=' ') ## end=" " here for i in range(n): print('{:13.4e}'.format(y[i]), end=' ') ## end=" " here print() m = Y.shape[0] try: n = Y.shape[1] except Exception: n = 1 if freq == 0: freq = m printHead(n) for i in range(0,m,freq): printLine(X[i], Y[i], n) if i != m - 1: printLine(X[m - 1], Y[m - 1], n) from numpy import array, zeros def F(x, y): F = zeros(2) F[0] = y[1] F[1] =(-y[1] - (y[0]/0.45) + 9.0)/2.0 return F x = 0.0 xStop = 10.0 H = 0.5 y = array([0.0, 0.0]) freq = 1 X, Y = bulStoer(F, x, y, xStop, H) printSoln(X, Y, freq)

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.