尝试在Python中实现Picard迭代

Question

我正在尝试写代码，用Python实现皮卡德迭代。特别是，我想在每一步迭代时绘制出近似值。这是我经过一天努力后得到的结果。

import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define the function f(u, t)
def f(u, t):
    return u

# Define the iteration function
def iterations(f, N, t_0, u_0, T, h):
    num_steps = int((T - t_0) / h)
    t = np.linspace(t_0, T, num_steps)

    # Compute the exact solution using odeint
    u_exact = odeint(f, u_0, t)
    plt.plot(t, u_exact, c = 'k', label="Exact solution")

    u = u_0*np.ones(num_steps)  # Initialize u
    
    # Perform iterations using the trapezoidal rule
    for j in range(0, N):
      plt.plot(t, u, label=f"Iteration {j}")
      integrand = lambda t: f(u, t)
      for i in range(num_steps):
          integral_result, _ = quad(integrand, t_0, t[i])
          u[i] = u_0 + integral_result


    plt.xlabel("t")
    plt.ylabel("u")
    plt.legend()
    plt.show()

iterations(f, 2, 0, 1, 2, 0.01)

我遇到的具体错误是只有长度为1的数组可以转换为Python标量：

类型错误 回溯（最近的调用最后） 在 <cell line: 35>() 33 plt.show() 34 ---> 35 iterations(f, 2, 0, 1, 2, 0.01)

2帧 /usr/local/lib/python3.10/dist-packages/scipy/integrate/_quadpack_py.py 在 _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points) 575 如果 points 是 None: 576 如果 infbounds == 0: --> 577 返回 _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit) 578 否则: 579 返回 _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

类型错误：只有长度为1的数组可以转换为Python标量

我想做的是在u[i]中记录与u(t_i)对应的值，而代码中的u实际上是数学中的第j次迭代。

Answer 1

完整的错误信息

TypeError                                 Traceback (most recent call last)
Cell In[1], line 35
     32     plt.legend()
     33     plt.show()
---> 35 iterations(f, 2, 0, 1, 2, 0.01)

Cell In[1], line 26, in iterations(f, N, t_0, u_0, T, h)
     24   integrand = lambda t: f(u, t)
     25   for i in range(num_steps):
---> 26       integral_result, _ = quad(integrand, t_0, t[i])
     27       u[i] = u_0 + integral_result
     30 plt.xlabel("t")

File ~\miniconda3\lib\site-packages\scipy\integrate\_quadpack_py.py:463, in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst, complex_func)
    460     return retval
    462 if weight is None:
--> 463     retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
    464                    points)
    465 else:
    466     if points is not None:

File ~\miniconda3\lib\site-packages\scipy\integrate\_quadpack_py.py:575, in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    573 if points is None:
    574     if infbounds == 0:
--> 575         return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
    576     else:
    577         return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

TypeError: only size-1 arrays can be converted to Python scalars

这个错误发生在调用 quad() 的时候，出现在几层之后。最后一层是一个编译过的函数，所以没有提供更多的细节。

但是如果我单独测试 integrand 这个函数

In [4]: f(np.ones(3),0)
Out[4]: array([1., 1., 1.])

它返回的数组大小和 u 参数是一样的。

如果我尝试把这个数组转换成 float 类型，就会出现你看到的错误信息

In [6]: float(f(np.ones(3),0))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[6], line 1
----> 1 float(f(np.ones(3),0))

TypeError: only size-1 arrays can be converted to Python scalars

我之前见过类似的错误，特别是在其他使用 quad 的情况下。不过 quad 的文档很明确，它只适用于标量值的函数。