Ошибка численного интегрирования с использованием Scipy.integrate - PullRequest
Новая
       97

Ошибка численного интегрирования с использованием Scipy.integrate

0 голосов
/ 02 мая 2020

Я пытаюсь выполнить следующую интеграцию. Код начинает выдавать ошибку, когда значение k превышает 4, а значение T c превышает 300. Сообщение об ошибке отображается внизу кода. Я провел несколько испытаний и обнаружил, что множители 2 и 0,5 в приведенных ниже утверждениях соответственно создают некоторую ошибку. Искал исправление для этой проблемы. 

B = (2 * (Ro ** 2 - ((Ro - a) ** 2))) ** 0.5
return [0, ((a - Ro) + ((Ro ** 2 - 0.5 * (z ** 2))) ** 0.5)]

import numpy as np
import scipy
from scipy import integrate
from scipy.integrate import nquad

Ro = 0.01
rat = 1
Ess = 207.8e9
k = 10
Tc = 900
Tm = 300

Ri = 0
a = rat * Ro

B = (2 * (Ro ** 2 - ((Ro - a) ** 2))) ** 0.5


def Cc55(y, z):
    R = (y ** 2 + z ** 2) ** 0.5
    Kc = 1.7 * (1 + 1.276 * (10 ** (-4)) * Tc + 6.648 * (10 ** (-8)) * (Tc ** 2))
    Km = 15.379 * (
        1
        - 1.264 * (10 ** (-3)) * Tm
        + 2.092 * (10 ** (-6)) * (Tm ** 2)
        - 7.223 * (10 ** (-10)) * (Tm ** 3)
    )
    a = (Kc - Km) / Km
    C = (
        1
        - (a / (k + 1))
        + (a ** 2 / (2 * k + 1))
        - (a ** 3 / (3 * k + 1))
        + (a ** 4 / (4 * k + 1))
        - (a ** 5 / (5 * k + 1))
    )
    d1 = (R - Ri) / (Ro - Ri)
    D = (
        d1
        - ((a * (d1 ** (k + 1))) / (k + 1))
        + (((a ** 2) * (d1 ** (2 * k + 1))) / ((2 * k) + 1))
        - (((a ** 3) * (d1 ** (3 * k + 1))) / ((3 * k) + 1))
        + (((a ** 4) * (d1 ** (4 * k + 1))) / ((4 * k) + 1))
        - (((a ** 5) * (d1 ** (5 * k + 1))) / ((5 * k) + 1))
    )
    b = D / C
    T = Tm + (Tc - Tm) * (b)
    EmT = (201.04 * (10) ** 9) * (
        1 + (3.079 * (10 ** (-4)) * T) - (6.534 * (10 ** (-7)) * (T ** 2))
    )

    EcT = (244.27 * (10) ** 9) * (
        1
        - (1.371 * (10 ** (-3)) * T)
        + (1.214 * (10 ** (-6)) * (T ** 2))
        - (3.681 * (10 ** (-10)) * (T ** 3))
    )
    VcT = (0.2882) * (1 + (1.133 * (10 ** (-4)) * T))

    VmT = (0.3262) * (
        1 - (2.002 * (10 ** (-4)) * T) + (3.797 * (10 ** (-7)) * (T ** 2))
    )
    Efg = (EcT - EmT) * ((d1) ** k) + EmT  # ((1.68e11-2.078e11)*((d1)**k) + 2.078e11)

    Vfg = (VcT - VmT) * ((d1) ** k) + VmT

    H = (4 * (Ro ** 2) - 2 * (z ** 2)) ** (0.5)

    L = (np.pi * y) / (2 * H)
    T = (y) / H
    M = (np.tan(L)) / L
    N1 = (1 - np.sin(L)) ** 4
    P = np.cos(L)
    F2 = ((M ** (0.5)) * (0.923 + 0.199 * (N1))) / P
    Q = Ro ** 2 - z ** 2
    S = (Ess) / (((Ro) ** 5) * Efg)

    return (16 * ((2) ** 0.5) * y * Q * (F2 ** 2)) * S


def bounds_z():
    return [0, B]


def bounds_y(z):
    return [0, ((a - Ro) + ((Ro ** 2 - 0.5 * (z ** 2))) ** 0.5)]


x4 = nquad(Cc55, [bounds_y, bounds_z])
print(x4[0])

C:\Users\arnab\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py:385: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  warnings.warn(msg, IntegrationWarning)
C:\Users\arnab\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py:385: IntegrationWarning: The integral is probably divergent, or slowly convergent.
  warnings.warn(msg, IntegrationWarning)
C:\Users\arnab\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py:385: IntegrationWarning: The algorithm does not converge.  Roundoff error is detected
  in the extrapolation table.  It is assumed that the requested tolerance
  cannot be achieved, and that the returned result (if full_output = 1) is 
  the best which can be obtained.
  warnings.warn(msg, IntegrationWarning)
...