平田智剛's questions

O que eu quero resolver

Resolvi um sistema de equações diferenciais ordinárias usando scipy.integrate.solve_ivpe comparei os resultados com minha implementação caseira de Runge-Kutta de quarta ordem (RK4).

Surpreendentemente, a precisão solve_ivp(usando RK45) é significativamente pior.

Alguém poderia me ajudar a entender por que isso pode estar acontecendo?

Descrição do problema

Simulei um movimento circular uniforme em um plano 2D com uma velocidade inicial de (10, 0) e uma aceleração centrípeta de magnitude 10. Teoricamente, esse sistema deveria descrever um círculo com raio de 10.

Plotei os resultados scipy.integrate.solve_ivp (method="RK45")em azul e aqueles da minha implementação caseira do RK4 em vermelho. O gráfico resultante é mostrado abaixo:

Como o RK4 tem precisão de 4ª ordem e o RK45 tem precisão de 5ª ordem, eu esperava que o RK45 tivesse um desempenho melhor e seguisse o círculo mais de perto. No entanto, os resultados do RK4 são muito superiores.

Código fonte relevante

from scipy.integrate import solve_ivp
import math
import matplotlib.pyplot as plt
import numpy as np

def get_a(t, r, v):
    # Accelerates perpendicular to the direction of motion with magnitude 10
    x, y = 0, 1
    
    direction_of_motion = math.atan2(v[y], v[x])
    direction_of_acceleration = direction_of_motion + math.pi / 2
    acceleration_magnitude = 10

    a_x = acceleration_magnitude * math.cos(direction_of_acceleration)
    a_y = acceleration_magnitude * math.sin(direction_of_acceleration)

    return np.array([a_x, a_y])

def get_v_and_a(t, r_and_v):
    # r_and_v == np.array([r_x, r_y, v_x, v_y])
    # returns    np.array([v_x, v_y, a_x, a_y])
    r = r_and_v[0:2]
    v = r_and_v[2:4]
    a = get_a(t, r, v)

    return np.concatenate([v, a])

# Simulation settings
initial_position = [0, 0]
initial_velocity = [10, 0]
end_time = 300
time_step = 0.01

# scipy.integrate.solve_ivp simulation (RK45)
initial_values = np.array(initial_position + initial_velocity)
time_points = np.arange(0, end_time + time_step, time_step)
result = solve_ivp(
    fun=get_v_and_a,
    t_span=[0, end_time],
    y0=initial_values,
    method="RK45",
    t_eval=time_points
)

scipy_ts = result.t
scipy_xs = result.y[0]
scipy_ys = result.y[1]

# Homemade RK4 simulation
my_ts, my_xs, my_ys = [], [], []
current_time = 0.0
step_count = 0
r = np.array(initial_position, dtype="float64")
v = np.array(initial_velocity, dtype="float64")
delta_t = time_step

while current_time <= end_time:
    my_ts.append(current_time)
    my_xs.append(r[0])
    my_ys.append(r[1])
    
    t = current_time + delta_t    
    r_and_v = np.concatenate([r, v])

    k0 = get_v_and_a(t -       delta_t, r_and_v                     )
    k1 = get_v_and_a(t - 1/2 * delta_t, r_and_v + 1/2 * k0 * delta_t)
    k2 = get_v_and_a(t - 1/2 * delta_t, r_and_v + 1/2 * k1 * delta_t)
    k3 = get_v_and_a(t                , r_and_v +       k2 * delta_t)

    step_count += 1
    current_time = step_count * delta_t
    r_and_v += (1/6 * k0 + 1/3 * k1 + 1/3 * k2 + 1/6 * k3) * delta_t
    r = r_and_v[0:2]
    v = r_and_v[2:4]

# Plot results
# 1. set up
size = max([abs(x) for x in scipy_xs] + [abs(y) for y in scipy_ys])*2.5
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_xlim(-size/2, size/2)
ax.set_ylim(-size/2, size/2)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_aspect("equal", adjustable="box")

# 2. plot
ax.plot(scipy_xs, scipy_ys, lw=0.1, color="blue",
    label="scipy.integrate.solve_ivp RK45")
ax.plot(my_xs, my_ys, lw=0.1, color="red", label="homemade code RK4")

# 3. time points
for i, t in enumerate(scipy_ts):
    if t % 5 == 0:
        ax.text(
            scipy_xs[i], scipy_ys[i], str(round(t)), fontsize=8,
            ha = "center", va = "center", color = "blue"
        )
    if t % 1 == 0 and (t <= 5 or 295 <= t): 
        ax.text(
            my_xs[i], my_ys[i], str(round(t)), fontsize=8,
            ha = "center", va = "center", color = "red"
        )

# 4. show
leg = ax.legend()
leg.get_lines()[0].set_linewidth(1)
leg.get_lines()[1].set_linewidth(1)
ax.grid(True)
plt.show()

Sobre o Código Caseiro

O que eu tentei

Li a documentação oficial , solve_ivpmas não consegui descobrir o que pode estar causando essa discrepância.

Estou trabalhando em um Python 3.13.1projeto usando mypy 1.14.0para verificação de tipo estático.
Tenho um módulo chamado module.pycom uma functionfunção que retorna um tipo com um nome muito longo,
Type_whose_name_is_so_long_that_we_do_not_want_to_call_it_over_and_over_again.
Para tornar o código mais legível, defini um alias de tipo Tno arquivo stub correspondente module.pyi.

Aqui está uma versão simplificada do meu código module.pyi:

T = Type_whose_name_is_so_long_that_we_do_not_want_to_call_it_over_and_over_again

def function()->T: pass

class Type_whose_name_is_so_long_that_we_do_not_want_to_call_it_over_and_over_again: pass

Quero evitar illegal_usage_of_T.pyque o seguinte use o Talias de tipo.

import module

foo:module.T = module.function()

O ideal é que, ao executar mypy illegal_usage_of_T.py, eu receba uma mensagem de erro indicando que o tipo Té indefinido.

O que eu tentei

Pesquisa Google

Pesquisei no Google por "mypy type alias only used in stub file", mas não consegui encontrar uma solução que impeça o mypy de reconhecer o type alias em outros módulos. Eu esperava que definir Tapenas no arquivo stub limitaria seu escopo, mas parece que o mypy consegue encontrar o type alias mesmo em outros módulos.

consertando código

Tentei várias abordagens, incluindo:

• Renomeando o alias do tipo: mudei Tpara _Tpara torná-lo menos provável de ser encontrado por outros módulos, mas isso não resolveu o problema.
• Usando if TYPE_CHECKING: Tentei definir condicionalmente o alias de tipo dentro de um if TYPE_CHECKINGbloco, mas isso também não impediu que o alias de tipo fosse usado em outros módulos.
• Limitando exportações: adicionei __all__ = ["function", "Type_whose_name_is_so_long_that_we_do_not_want_to_call_it_over_and_over_again"]ao module.pyiarquivo para controlar explicitamente quais nomes são exportados, mas o alias do tipo Tainda estava acessível.

Aqui está o código modificado para module.pyi:

from typing import TYPE_CHECKING

if TYPE_CHECKING:
    _T = Type_whose_name_is_so_long_that_we_do_not_want_to_call_it_over_and_over_again

def function()->_T: pass

class Type_whose_name_is_so_long_that_we_do_not_want_to_call_it_over_and_over_again: pass

__all__ = ["function", "Type_whose_name_is_so_long_that_we_do_not_want_to_call_it_over_and_over_again"]

E aqui está o illegal_usage_of_T.pyarquivo:

import module

foo:module._T = module.function()  # Still works

Eu esperava que pelo menos uma dessas abordagens impedisse illegal_usage_of_T.pyo acesso ao _Talias de tipo, mas nenhuma delas funcionou.

Embora a assinatura do método em Subseja compatível com Super, mypy rejeita a substituição: Signature of "method" incompatible with supertype "Super".

estou usando

• python 3.13.1
• meupy 1.14.0

Primeiro, fiz o seguinte test.pyi.

from typing import overload

class Super:
    def method(self, arg:Other|Super)->Super: pass

class Sub(Super):
    @overload
    def method(self, arg:Other|Sub)->Sub: pass
    @overload
    def method(self, arg:Super)->Super: pass

class Other: pass

Então, quando executei mypy test.pyina linha de comando, o mypy produziu o seguinte diagnóstico:

test.pyi:7: error: Signature of "method" incompatible with supertype "Super"  [override]
test.pyi:7: note:      Superclass:
test.pyi:7: note:          def method(self, arg: Other | Super) -> Super
test.pyi:7: note:      Subclass:
test.pyi:7: note:          @overload
test.pyi:7: note:          def method(self, arg: Other | Sub) -> Sub
test.pyi:7: note:          @overload
test.pyi:7: note:          def method(self, arg: Super) -> Super
Found 1 error in 1 file (checked 1 source file)

Verifiquei o tipo de E/S de ambos Super.methode Sub.methode descobri que não há nenhum padrão que viole o LSP (Princípio de Substituição de Liskov).

Sub.methodLata sobrecarregada

• entrada argdo tipo Other|Super(= Other|Sub+ Super) e
• tipo de saída Super(= Sub+ Super).

O tipo de entrada e saída acima corresponde à assinatura de Super.method.

Então, não tenho ideia do que ser Signature of "method" incompatible with supertype "Super".

A seguir está a tabela de E/S de method.

*O LSP exige que o tipo de retorno de um submétodo seja mais estreito ou igual ao tipo de retorno do supermétodo.

** Sub.methodretorna

• Subquando Othera entrada e
• Superquando Supera entrada,

então ele retorna Sub|Superquando Other|Supera entrada.
Sub|Supersignifica Super.

Como você pode ver na tabela acima, não há padrões que violem o LSP.

Então, acho que a mensagem de erro mypy Signature of "method" incompatible with supertype "Super"está incorreta.

Meu código está errado?

Além disso, se meu código não estiver errado e a mensagem de erro do mypy estiver errada, onde posso perguntar?

PS Uma maneira rápida de esconder o erro.

Embora não seja uma solução completa, encontrei uma maneira simples de esconder o erro.

from typing import overload

class Super:
    def method(self, arg:Other|Super)->Super: pass

class Sub(Super):
    @overload
    def method(self, arg:Other|Sub)->Sub: pass
    @overload
    def method(self, arg:Super)->Super: pass
    # Start of hiding error
    @overload
    def method( # type: ignore[overload-cannot-match]
               self, arg:Other|Super)->Super: pass
    # End of hiding error

class Other: pass

Como uma última sobrecarga, adicionei um Sub.methodcom a mesma assinatura exata de Super.method. No entanto, quando esse problema for resolvido em uma versão futura, isso significará que o código que adicionei não será alcançado e devemos obter um [overload-cannot-match]erro. Portanto, adicionei # type: ignore[overload-cannot-match]para ignorar esse erro com antecedência.

(À primeira vista, pode parecer que estou simplesmente silenciando erros com type: ignore, mas isso não é relevante para o mypy no momento. Isso é apenas um impedimento contra erros futuros.)