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

Question

在下面的代码中，为什么我一直收到这个错误：

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)

Answer 1

1. “def integration_BulirschStoeir(F, x, y, xStop, tol)”中存在缩进问题（从“kMax = 51”到函数结束）。因此，在 "X, Y = bulStoer(F, x, y, xStop, H)" 之后 Y 中没有新数据，并且 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)