Метод Эйлера для разных размеров шагов. Как изменить код алгоритма для учета различных значений размера шага?

Question

У меня есть алгоритм для задачи численного дифференциального уравнения, называемый методом Эйлера. По сути, метод Эйлера приближает решение дифференциальных уравнений. Моя функция работает для одного размера шага (значение h), но я пытаюсь изменить код, чтобы разрешить мне l oop для 3 различных значений h (путем изменения h с одного значения на список возможные значения). Тем не менее, функция, которую я написал, не зацикливается на моих значениях. Я новичок в python и ранее использовал R. Может кто-нибудь показать мне, как это сделать правильно.

Мой код, который работает для одного значения шага h:

from math import exp # exponential function

dy = lambda x,y: x*y
f = lambda x: exp(x**2/2) # analytical solution function

x = 0 # Intial value X_0
xn = 2 # Final Value
y = 1 # value of y(x0)
h = 0.2 # stepsize
n = int((xn-x)/h)

print ('x \t\t y (Euler h={}) \t y (analytical)'.format(h))
print ('%f \t %f \t %f'% (x,y,f(x)))
for i in range(n):
    y += dy(x, y)*h
    x += h
    print ('%f \t %f \t %f'% (x,y,f(x)))


x        y (Euler h=0.5) y (analytical)
0.000000     1.000000    1.000000
0.500000     1.000000    1.133148
1.000000     1.250000    1.648721
1.500000     1.875000    3.080217
2.000000     3.281250    7.389056

Я хотел бы изменить h на h=[0.01,0.2,0.5] и иметь значения, чтобы затем создать график, показывающий аналитическое решение и решения метода Эйлера при различных значениях размера шага.

Еще раз прошу прощения, если это простой вопрос. Я новичок в программировании в python и продолжаю делать некоторые ошибки, ниже моя лучшая попытка на данный момент. Я еще не сохранил свои значения x в контейнере, так как моя функция не зацикливалась на значениях h. Я пытаюсь написать вложенное значение для l oop, где внешний l oop перебирает значения h, сохраняет значения и отображает их в виде строки, затем переходит ко второму значению h и делает то же самое на В конце значения могут быть размещены на одном графике.

# Improved to allow plotting different values
import matplotlib.pyplot as plt 
import numpy as np
from math import exp # exponential function

dy = lambda x,y: x*y
f = lambda x: exp(x**2/2) # analytical solution function
x = 0
xn = 2
y = 1
# Container for step sizes
h = [0.5,0.2,0.1]

# Container to store the x values at each stepsize
# X =np.zeros((3,))

print ('x \t\t y (Euler) \t y (analytical)')
print ('%f \t %f \t %f'% (x,y,f(x)))
for j in range(0,len(h),1):
    n = int((xn-x)/h[j])
    for i in range(n):
        y += dy(x, y)*h[j]
        x += h[j]
        print ('%f \t %f \t %f'% (x,y,f(x)))
    plt.plot(x,y)

plt.show()


x        y (Euler)   y (analytical)
0.000000     1.000000    1.000000
0.500000     1.000000    1.133148
1.000000     1.250000    1.648721
1.500000     1.875000    3.080217
2.000000     3.281250    7.389056

Так что вопрос на самом деле пытается создать метод Эйлера для разных размеров шагов т.е. "Как можно изменить нашу функцию на l oop по списку и отобразить результаты с помощью matplotlib"?

Answer 1

Вы допустили небольшую ошибку, и вам нужно сохранить результаты в контейнере, если вы хотите построить их. Я немного переписал твой код. Сначала я дам вам полный код, прежде чем обсуждать, что не так с вашим кодом. Может быть, вы сами заметите ошибки. Я также добавил расчет аналитического решения и некоторые другие небольшие улучшения, которые должны вам понравиться. Итак, вот код:

import matplotlib.pyplot as plt
import numpy as np
from math import exp  # exponential function

dy = lambda x, y: x * y
f = lambda x: exp(x ** 2 / 2)  # analytical solution function
x_final = 2

# analytical solution
x_a = np.arange(0, x_final, 0.01)
y_a = np.zeros(len(x_a))
for i in range(len(x_a)):
    y_a[i] = f(x_a[i])
plt.plot(x_a, y_a, label="analytical")

# Container for step sizes
h = [0.5, 0.2, 0.1]


for j in range(len(h)):
    x = 0
    y = 1
    print("h = " + str(h[j]))
    print("x \t\t y (Euler) \t y (analytical)")
    print("%f \t %f \t %f" % (x, y, f(x)))

    n = int((x_final - x) / h[j])

    x_n = np.zeros(n + 1)
    y_n = np.zeros(n + 1)
    x_n[0] = x
    y_n[0] = y

    for i in range(n):
        y += dy(x, y) * h[j]
        x += h[j]
        print("%f \t %f \t %f" % (x, y, f(x)))
        x_n[i + 1] = x
        y_n[i + 1] = y

    plt.plot(x_n, y_n, "x-", label="h=" + str(h[j]))


plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

На моем компьютере показано следующее:

Обратите внимание, что я переименовал вашу переменную xn в x_final чтобы избежать путаницы с введенными мною переменными. Как указано выше, вам необходимо хранить каждое из ваших значений x и y в контейнере. Для этого я использовал NumPy массивы, но вы также можете использовать список. Этот

n = int((x_final - x) / h[j])

x_n = np.zeros(n + 1)
y_n = np.zeros(n + 1)
x_n[0] = x
y_n[0] = y

просто создает 2 массива нулей, размер которых равен числу подэтапов +1. Затем я устанавливаю первые значения равными начальным значениям. Это должно быть внутри l oop сверх h, так как число подшагов n отличается для каждого h.

В конце вашего i -l oop, я просто пишу текущие значения x и y для правильной позиции в массивах.

for i in range(n):
    y += dy(x, y) * h[j]
    x += h[j]
    print("%f \t %f \t %f" % (x, y, f(x)))
    x_n[i + 1] = x
    y_n[i + 1] = y

Вместо вызова plt.plot с x и y, который просто отображает одну точку, так как они являются sclars, вам нужно передать массивы в функцию:

plt.plot(x_n, y_n, "x-", label="h=" + str(h[j]))

Я также добавил метку, которая будет отображаться в легенде, и изменил тип линии на "x-".

Вы допустили одну ошибку, из-за которой ваши i l oop были выполнены только для первых h, заключались в том, что вы не сбросили x и y к их начальным значениям. Таким образом, ваш n всегда был 0 после первого запуска внешнего l oop.

Конечно, есть несколько вещей, которые вы можете оптимизировать, например, использовать что-то вроде

for h in h_list:
   ...

который был бы немного более читабельным, чем всегда, используя h[j] вместо просто h, но я думаю, что пока этого достаточно. ;)