Исследование пересечения искусственного интеллекта и физики: роль ChatGPT в генерации кода 2023

Представьте себе мир, где машины могут генерировать код для решения сложных проблем в окружающем нас физическом мире. ChatGPT,  Natural Language Processor (NLP), который записывает человекоподобные ответы из пользовательских подсказок ввода, может сделать именно это.

В этой статье я собираюсь показать вам, как это делается. Прямо сейчас любой желающий может воспользоваться исследовательской версией ChatGPT — вам просто нужно зайти на веб-сайт OpenAIs и зарегистрировать учётную запись, чтобы опробовать данный инструмент.

Многое происходит под капотом ChatGPT, и я не собираюсь пытаться объяснять вам это здесь (OpenAI даёт подробный обзор того, как работает технология на своём веб-сайте).

Вместо этого я собираюсь проработать серию примеров решения задач, основанных на физике, каждая из которых имеет возрастающую сложность. К концу у вас должно быть хорошее представление о сильных и слабых сторонах ChatGPT применительно к такого рода проблемам.

Мотивация

Прежде чем углубляться в примеры, неплохо бы начать с вопроса, почему мы должны быть заинтересованы в использовании искусственного интеллекта для написания подобного кода. Будучи человеком, который большую часть своей карьеры проработал в индустрии программного обеспечения для автоматизированного проектирования (CAE), я считаю, что это может стать масштабным прорывным приложением для новых технологий искусственного интеллекта.

Согласно Fortune Business Insights, объём мирового рынка CAE-инжиниринга, по прогнозам, вырастет с 6,25 млрд долларов США до 12,16 млрд долларов США к 2029 году. CAE используется в жизненном цикле разработки стольких продуктов и процессов, что трудно переоценить ценность, которую программные продукты CAE предоставляют, часто косвенно, стольким людям по всему миру каждый божий день. В любом месте продукта или процесса, где существует какая-либо структурная, динамическая, текучая, тепловая или электромагнитная среда, которую необходимо понять или оптимизировать, при разработке, скорее всего, использовалась та или иная форма CAE. Самолёты, автомобили и персональные электронные устройства – всё это продукты, в разработке которых применяется CAE. Переработка добытых руд, литьё пластмасс под давлением и переработка вторичных отходов – всё это процессы, которые были оптимизированы с помощью CAE. Инновации и усовершенствования в программных продуктах CAE многократно повышают производительность людей, использующих эти инструменты. Если у вас на выполнение моделирования уходит вдвое меньше времени, то инженеры, использующие это программное обеспечение, могут опробовать вдвое больше проектов за тот же промежуток времени. Если вы можете ввести какую-либо форму новой физики или парного моделирования (например, взаимодействие жидкости и структуры), тогда инженеры смогут учитывать мультифизику в своей среде моделирования, что ещё больше упростит процесс разработки.

Программные продукты CAE сложны и часто адаптируются к одному конкретному приложению — анализ конечных элементов (FEA) используется для анализа структур, а вычислительная гидродинамика (CFD) используется, например, для моделирования потоков жидкости. Одна вещь, которая у всех них является общей – это физическая модель. Эти физические модели обычно состоят из некоторого набора управляющих уравнений и эффективного, с вычислительной точки зрения, алгоритма их решения. Именно разработку этих физических моделей с помощью ChatGPT я исследую здесь.

Я попрошу ChatGPT написать код для моделирования следующих физических систем:

1. Двойной маятник
2. Теплопередача
3. Уравнения Навье – Стокса

Я буду разбивать каждый из этих примеров на объяснение, подсказку, которую я использую, ответ ChatGPT и первые впечатления.

Двойной маятник

Двойной маятник – это высокодинамичная система, поэтому его моделирование должно дать несколько приятных результатов, изменяющихся во времени. Эта система состоит из двух конечностей, закреплённых на одном конце. Суставы, соединяющие конечности, не подвержены трению. При некоторых начальных условиях положение каждой ветви, образующей двойной маятник, может быть перенесено вперёд во времени. Двойной маятник очень чувствителен к этим начальным условиям. Мы можем думать о двойном маятнике как о простой системе многотелесной динамики (MBD). Полноценные MBD-коды могут решать гораздо более сложные задачи, такие как движение сельскохозяйственной техники или динамика внедорожника на сложной местности. Давайте посмотрим, как ChatGPT справляется с моей подсказкой для решения этой проблемы.

Подсказка

“Напишите код на Python, который имитирует двойной маятник. Постройте график и анимируйте результаты с помощью matplotlib.”

Ответ

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# Physical constants
g = 9.8  # acceleration due to gravity
l1 = 1  # length of pendulum 1
l2 = 1  # length of pendulum 2
m1 = 1  # mass of pendulum 1
m2 = 1  # mass of pendulum 2

# Initial conditions
theta1_0 = np.pi / 2
theta2_0 = np.pi / 2
omega1_0 = 0
omega2_0 = 0

# Time step and number of steps
dt = 0.01
steps = 2000

# Calculate the derivatives of theta1 and theta2
def derivatives(state, t):
    theta1, theta2, omega1, omega2 = state
    dydx = np.zeros_like(state)
    dydx[0] = omega1
    dydx[1] = omega2

    delta = theta2 - theta1
    den1 = (m1 + m2) * l1 - m2 * l1 * np.cos(delta) * np.cos(delta)
    dydx[2] = (m2 * l1 * omega1 * omega1 * np.sin(delta) * np.cos(delta) +
               m2 * g * np.sin(theta2) * np.cos(delta) +
               m2 * l2 * omega2 * omega2 * np.sin(delta) -
               (m1 + m2) * g * np.sin(theta1)) / den1

    dydx[3] = (-m2 * l2 * omega2 * omega2 * np.sin(delta) * np.cos(delta) +
               (m1 + m2) * g * np.sin(theta1) * np.cos(delta) -
               (m1 + m2) * l1 * omega1 * omega1 * np.sin(delta) -
               (m1 + m2) * g * np.sin(theta2)) / (l2 * den1)
    return dydx

# Integrate the equations of motion using the RK4 method
def integrate(state, t):
    k1 = dt * derivatives(state, t)
    k2 = dt * derivatives(state + 0.5 * k1, t + 0.5 * dt)
    k3 = dt * derivatives(state + 0.5 * k2, t + 0.5 * dt)
    k4 = dt * derivatives(state + k3, t + dt)
    return state + (k1 + 2 * k2 + 2 * k3 + k4) / 6

# Calculate the positions of the pendulums
def calculate_positions(state):
    x1 = l1 * np.sin(state[0])
    y1 = -l1 * np.cos(state[0])
    x2 = x1 + l2 * np.sin(state[1])
    y2 = y1 - l2 * np.cos(state[1])
    return x1, y1, x2, y2

# Set up the initial state
state = np.array([theta1_0, theta2_0, omega1_0, omega2_0])
t = np.arange(0, steps * dt, dt)
x1_list = []
y1_list = []
x2_list = []
y2_list = []

for i in range(steps):
    x1, y1, x2, y2 = calculate_positions(state)
    x1_list.append(x1)
    y1_list.append(y1)
    x2_list.append(x2)
    y2_list.append(y2)
    state = integrate(state, t[i])

# Update the animation for each frame
def animate(i):
    x1 = [0, x1_list[i]]
    y1 = [0, y1_list[i]]
    x2 = [x1_list[i], x2_list[i]]
    y2 = [y1_list[i], y2_list[i]]
    x3 = x1_list[:i]
    y3 = y1_list[:i]
    x4 = x2_list[:i]
    y4 = y2_list[:i]
    line1.set_data(x1, y1)
    line2.set_data(x2, y2)
    line3.set_data(x3, y3)
    line4.set_data(x4, y4)
    return line1, line2, line3, line4

# Set up the figure for the animation
fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
ax.set_aspect('equal')
line1, = ax.plot([], [], '-', lw=2, color='red')
line2, = ax.plot([], [], '-', lw=2, color='blue')
line3, = ax.plot([], [], '-', lw=1, color='red')
line4, = ax.plot([], [], '-', lw=1, color='blue')

# Initialize the animation
def init():
    line1.set_data([], [])
    line2.set_data([], [])
    line3.set_data([], [])
    line4.set_data([], [])
    return line1, line2, line3, line4

# Run the animation
ani = animation.FuncAnimation(fig, animate, np.arange(1, len(y1_list)),
                              interval=25, blit=True, init_func=init)

plt.show()

Первые впечатления

ChatGPT дал очень разумный, хорошо прокомментированный ответ. Система решается способом пространства состояний, где вычисляются производные состояния, а затем используются для обновления будущего состояния системы. Однако, что впечатляет меня больше всего, так это функция интеграции. Если бы я сам впервые попытался решить эту проблему, я бы, вероятно, предпочёл использовать метод интегрирования Эйлера, но нет, ChatGPT сразу перешёл к использованию метода Рунге-Кутты (RK4), который обычно даёт меньшие ошибки за счёт более дорогостоящих вычислений. Он сделал именно то, что я просил, и даже больше (используя RK4, когда я явно не запрашивал его об этом). Отличное начало. Почему бы нам не попробовать что-нибудь более сложное?

Теплопередача

Мы видели, как ChatGPT успешно произвёл симуляцию системы, изменяющейся во времени. Теперь давайте посмотрим, как он справится с системой, изменяющейся в пространстве и во времени. В этом примере я попрошу ChatGPT произвести имитацию одномерной теплопередачи вдоль стальной балки с фиксированной температурой на каждом конце. Решение этой проблемы требует решения уравнения теплопроводности, уравнения в частных производных (PDE), которое является неотъемлемой частью множества различных типов программного обеспечения CAE. Теплопередача рассматривается во всех видах различных инженерных приложений. Способность решать такого рода проблемы важна, поэтому давайте посмотрим, как ChatGPT справится с этим.

Подсказка

“На Python напишите код для моделирования одномерной теплопередачи вдоль стальной балки длиной 0,1 м и площадью 0,01 м2, используя уравнение теплопроводности. Один конец балки должен быть зафиксирован на 373,15 Кельвина, а другой – на 273,15 Кельвина. Чтобы визуализировать результаты, постройте и анимируйте с течением времени линейный график, используя FuncAnimation в matplotlib. Включите временную метку на графике.”

Ответ

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# Constants
L = 0.1    # length of beam in meters
A = 0.01   # area of beam in square meters
k = 60     # thermal conductivity of steel in W/(m*K)
rho = 7850 # density of steel in kg/m^3
Cp = 480   # specific heat capacity of steel in J/(kg*K)
T0 = 373.15 # temperature at one end of the beam in Kelvin
T1 = 273.15 # temperature at the other end of the beam in Kelvin
dx = 0.002 # distance step size in meters
dt = 0.001  # time step size in seconds
t_end = 2 # end time in seconds

# Derived constants
alpha = k / (rho * Cp * A)
N = int(L/dx)
M = int(t_end/dt)

# Initial conditions
T = np.zeros((N, M))
T[:,0] = T1

# Boundary conditions
T[0,:] = T0
T[-1,:] = T1

# Heat equation solver
for i in range(0, M-1):
    for j in range(1, N-1):
        T[j,i+1] = T[j,i] + alpha * (T[j+1,i] - 2*T[j,i] + T[j-1,i]) * dt / dx**2

# Plotting function
fig, ax = plt.subplots()
line, = ax.plot([], [], lw=2)
ax.set_xlim(0, L)
ax.set_ylim(273, 373)
ax.set_xlabel('Distance (m)')
ax.set_ylabel('Temperature (K)')
ax.set_title('1D Heat Transfer along a Steel Beam')
time_text = ax.text(0.05, 0.95, '', transform=ax.transAxes)

# Animation function
def animate(i):
    line.set_data(np.linspace(0, L, N), T[:,i])
    time_text.set_text('Time = {} s'.format(round(i*dt, 3)))
    return line, time_text

# Create animation
anim = FuncAnimation(fig, animate, frames=M, interval=1, blit=True)

# Show plot
plt.show()

Первые впечатления

Дисклеймер: Этот ответ – не совсем то, что дал мне ChatGPT. Исходный ответ установил dx = 0.001 и dt = 0.01, что привело к сбою моделирования. Это происходит, если сетка слишком мелкая по сравнению с размером временного шага, поэтому я скорректировал dx и dt, пока не получил стабильное моделирование.

В подсказке для этого случая я уточнил длину и площадь поперечного сечения балки, но что касается свойств материала, я просто сказал, что она сделана из стали. Это было сделано специально, мне было интересно посмотреть, как ChatGPT справится с этим, и он проделал хорошую работу по настройке свойств материала. Значения теплопроводности, плотности и удельной теплоёмкости находятся в диапазоне, который можно было бы ожидать от стального сплава.

Обращая внимание на алгоритм уравнения теплопроводности, это было сделано во вложенном цикле for, поэтому это не самый эффективный с точки зрения вычислений метод, но он решает уравнение и даёт ожидаемый результат — со временем тепло рассеивается вдоль луча и в конечном итоге приводит к равновесному состоянию с постоянным температурным градиентом вдоль луча. За исключением проблемы зависимости размера сетки от временного шага и, возможно, эффективности решения уравнения теплопроводности, ChatGPT снова дал хороший отклик.

Уравнения Навье – Стокса

Пока всё хорошо. Я поставил две сложные задачи ChatGPT, и оба раза он дал мне хорошие ответы. В последнем примере я собираюсь посмотреть, смогу ли я заставить ChatGPT написать полноценный CFD-решатель для уравнений Навье-Стокса. Я заинтересован в получении полного развёрнутого решения.

Если вы можете написать код для решения этой задачи, значит, вы хорошо разбираетесь в уравнениях Навье-Стокса и численных методах, необходимых для их решения. Если вы заинтересованы в более глубоком понимании этого материала, ознакомьтесь с удивительной книгой Лорены Барба “12 Steps to Navier-Stokes“, в которой собраны части, необходимые для написания конечно-разностного CFD-решателя. CFD используется в различных типах применений, где необходимо понимать поведение жидкости (самолёты, дорожные транспортные средства, горение и так далее).

Это также особенно трудная область для освоения из-за сложного поведения текучих сред и численных методов, необходимых для решения этой категории задач.

Подсказка

“На Python напишите конечно-разностный вычислительный решатель гидродинамики и решите уравнение Навье-Стокса. Отобразите результаты с помощью matplotlib, используя “pcolor”.

Ответ

import numpy as np
import matplotlib.pyplot as plt

# Parameters
Lx, Ly = 1.0, 1.0  # Domain size
nx, ny = 101, 101  # Number of grid points
dx, dy = Lx/(nx-1), Ly/(ny-1)  # Grid spacing
Re = 100.0  # Reynolds number
dt = 0.001  # Time step
T = 10.0  # Final time

# Grid
x = np.linspace(0.0, Lx, nx)
y = np.linspace(0.0, Ly, ny)
X, Y = np.meshgrid(x, y)

# Initial conditions
u = np.zeros((ny, nx))
v = np.zeros((ny, nx))
p = np.zeros((ny, nx))

# Boundary conditions
u[:, 0] = 0.0  # Left wall
u[:, -1] = 0.0  # Right wall
u[0, :] = 0.0  # Bottom wall
u[-1, :] = 1.0  # Top wall
v[:, 0] = 0.0  # Left wall
v[:, -1] = 0.0  # Right wall
v[0, :] = 0.0  # Bottom wall
v[-1, :] = 0.0  # Top wall

# Solver
for n in range(int(T/dt)):
    # Pressure Poisson equation
    b = np.zeros((ny, nx))
    b[1:-1, 1:-1] = (1.0/dt)*((u[1:-1, 2:]-u[1:-1, 0:-2])/(2.0*dx) + 
                               (v[2:, 1:-1]-v[0:-2, 1:-1])/(2.0*dy))
    for q in range(100):
        pn = p.copy()
        p[1:-1, 1:-1] = ((pn[1:-1, 2:]+pn[1:-1, 0:-2])*dy**2 +
                         (pn[2:, 1:-1]+pn[0:-2, 1:-1])*dx**2 -
                         b[1:-1, 1:-1]*dx**2*dy**2)/(2.0*(dx**2+dy**2))
        p[:, -1] = p[:, -2]  # Right wall
        p[0, :] = p[1, :]  # Bottom wall
        p[:, 0] = p[:, 1]  # Left wall
        p[-1, :] = p[-2, :]  # Top wall

    # Velocity update
    u[1:-1, 1:-1] = (u[1:-1, 1:-1] -
                     u[1:-1, 1:-1]*dt/dx*(u[1:-1, 1:-1]-u[1:-1, 0:-2]) -
                     v[1:-1, 1:-1]*dt/dy*(u[1:-1, 1:-1]-u[0:-2, 1:-1]) -
                     dt/(2.0*Re*dx)*(p[1:-1, 2:]-p[1:-1, 0:-2]) +
                     dt/Re*(u[1:-1, 2:]-2.0*u[1:-1, 1:-1]+u[1:-1, 0:-2])/dx**2 +
                     dt/Re*(u[2:, 1:-1]-2.0*u[1:-1, 1:-1]+u[0:-2, 1:-1])/dy**2)
    v[1:-1, 1:-1] = (v[1:-1, 1:-1] -
                     u[1:-1, 1:-1]*dt/dx*(v[1:-1, 1:-1]-v[1:-1, 0:-2]) -
                     v[1:-1, 1:-1]*dt/dy*(v[1:-1, 1:-1]-v[0:-2, 1:-1]) -
                     dt/(2.0*Re*dy)*(p[2:, 1:-1]-p[0:-2, 1:-1]) +
                     dt/Re*(v[1:-1, 2:]-2.0*v[1:-1, 1:-1]+v[1:-1, 0:-2])/dx**2 +
                     dt/Re*(v[2:, 1:-1]-2.0*v[1:-1, 1:-1]+v[0:-2, 1:-1])/dy**2)

    # Boundary conditions
    u[:, 0] = 0.0  # Left wall
    u[:, -1] = 0.0  # Right wall
    u[0, :] = 0.0  # Bottom wall
    u[-1, :] = 1.0  # Top wall
    v[:, 0] = 0.0  # Left wall
    v[:, -1] = 0.0  # Right wall
    v[0, :] = 0.0  # Bottom wall
    v[-1, :] = 0.0  # Top wall

# Plotting
fig = plt.figure(figsize=(8, 6))
plt.pcolor(X, Y, np.sqrt(u**2 + v**2), cmap='GnBu_r')
colorbar = plt.colorbar()
colorbar.set_label('u [ms$^{-1}$]')
skip = (slice(None, None, 4), slice(None, None, 4))
plt.quiver(X[skip], Y[skip], u[skip], v[skip], color='white')
plt.title('2D Lid-Driven Cavity Flow')
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.show()

Первые впечатления

Дисклеймер: Мне действительно пришлось изменить пару вещей из ответа. Первая из них заключалось в том, что ChatGPT установил для крышки граничное условие нулевого давления. Обычно это условие выхода, и не имеет смысла включать его в процесс решения данной задачи.

Я исправил это, чтобы граничное условие нулевого градиента давления применялось ко всем стенам.

Я увеличил время моделирования до 10 секунд, чтобы дать больше времени для разработки стационарного решения, и добавил график дрожания, чтобы помочь визуализировать поток.

Говоря о визуализации, я также изменил цветовую карту с jet на GnBu_r. Существует множество причин, по которым jet является плохим выбором для визуализации данных. Натаниэль Смит выступил с замечательным докладом именно на эту тему на SciPy2015 https://www.youtube.com/watch?v=xAoljeRJ3lU, который я настоятельно рекомендую посмотреть.

Глядя на код, в нём есть всё, что вы ожидали бы увидеть для решения такого рода задач с использованием метода конечных разностей. Он определил сетку для решения задачи, константы, такие как число Рейнольдса, и определил граничные условия.

Уравнение обеспечивает способ вычисления поля давления из поля скорости и гарантирует, что поле скорости удовлетворяет сохранению массы и импульса. Это делается во внутреннем цикле. Затем поле скоростей обновляется, и в конце каждого шага применяются граничные условия. ChatGPT также вывел результаты в простой, читаемой форме (которую я впоследствии изменил по своему вкусу).

В отличие от примера с теплопередачей, этот код векторизован, что хорошо сказывается на производительности. Таким образом, из двухстрочного запроса ChatGPT выдал около 80 строк кода на Python, который может решать уравнения Навье-Стокса!

Заключение

В этой статье я подробно рассмотрел, как ChatGPT может быть использован для помощи в написании кода для решения реальных инженерных задач, кульминацией которого является создание CFD-решателя на основе конечных разностей для решения уравнения Навье-Стокса. ChatGPT продемонстрировал отличные способности в понимании подсказок, которые я ему давал.

Он далёк от совершенства (смотрите, сколько раз мне приходилось вмешиваться, когда он делал что-то странное), и это необходимо подчеркнуть. Я попросил ChatGPT решить проблемы в домене, с которым я знаком, и поэтому, когда он выдал эти странные ответы, я знал, что я ищу и что мне нужно изменить, чтобы устранить проблемы. Это предупреждение о том, что при использовании такого рода инструментов следует действовать очень осторожно. Они могут “укусить” вас, а вы можете даже не осознавать этого.

В начале этой статьи я утверждал, что новые технологии искусственного интеллекта, такие как ChatGPT, могут быть чрезвычайно разрушительными, когда используются для оказания помощи инженерам в решении сложных реальных проблем.

Хотя примеры, показанные в этой статье, далеко не так сложны, как производственный CAE-код, я надеюсь, что они послужат демонстрацией того, что возможно с помощью ChatGPT для решения задач, основанных на физике. Нам ещё так много предстоит узнать о том, как можно использовать эти инструменты.

Например, я уже начал использовать ChatGPT для помощи в написании модульных тестов для pySAMSS, проекта с открытым исходным кодом.

Куда эта технология пойдёт дальше, можно только догадываться. Если бы вы сказали мне 5 лет назад, что появится такая штука, которая сможет писать код для решения сложных физических процессов, я бы вам не поверил. Мне, например, не терпится увидеть, куда денутся эти инструменты в ближайшие 5 лет!