尝试在Python中实现Picard迭代

Question

我正在尝试编写代码来在 Python 中实现 Picard 的迭代。特别是，我想绘制迭代每一步的近似值。这是我经过一整天的努力得到的结果。

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)

我在这里遇到的具体错误是

only length-1 arrays can be converted to Python scalars

。

我想做的是在

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

文档很清楚它只集成了标量值函数。