Краевая задача метод конечных разностей python

Решение сингулярной двухточечной краевой задачи

Нужно найти y(x) методом конечных разностей (изобразить график) для . p, q и f заданы.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

import numpy as np import matplotlib.pyplot as plt # сетка узлов x_i, p(x_i)=p_i, q(x_i)=x_i, f(x_i)=f_i def find_params(a, b, h): x = np.arange(a, b + h, h) p = x[1:-1] * np.sin(x[1:-1]) + x[1:-1] ** 2 q = np.sqrt(x[1:-1] + x[1:-1] ** 2) f = np.sqrt(x[1:-1] ** 2 - x[1:-1] ** 3) + x[1:-1] return x, p, q, f a = 0 b = 1 n = 80 eps_list = [1, 0.1, 0.01, 0.001] a1 = 1 a2 = 0 b1 = 1 b2 = 0 g1 = 6 g2 = 9 for i, eps in enumerate(eps_list): fig, ax = define_plot(y_phi=0) # задаёт свойства графика h = (b - a) / n x, p, q, f = find_params(a, b, h) diag_dn = np.hstack((eps - p * h / 2, -b1)) diag_0 = np.hstack((a2 * h +a1, q * (h ** 2) - 2 * eps, b2 * h + b1)) diag_up = np.hstack((-a1, eps + p * h / 2)) W = np.diag(diag_dn, k=-1) + np.diag(diag_0, k=0) + np.diag(diag_up, k=1) F = np.hstack((g1 * h, -f * (h ** 2), g2 * h)) Y = tridiagonal_algorithm(W, F) # решение СЛАУ WY=F ax.plot(x, Y, '-o', color=colors_list[i], label='$\\epsilon=$' + str(eps)) ax.legend(loc='lower right') fig.savefig(str(i)+'_j', dpi=300)

Проблема:
Для всех эпсилон, кроме последнего, всё работает хорошо, но для творится что-то, что я не вполне понимаю. При разных n функция то возрастающая, то убывающая. Например, для получаю убывающую функцию, то же для (кстати, величина максимума на графике тоже меняется), а для функция уже возрастает. Кто-нибудь может подсказать, с чем это связано?

Источник

Saved searches

Use saved searches to filter your results more quickly

You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session. You switched accounts on another tab or window. Reload to refresh your session.

KrivosheyevaMarina/SummerPractic

This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

Sign In Required

Please sign in to use Codespaces.

Launching GitHub Desktop

If nothing happens, download GitHub Desktop and try again.

Launching GitHub Desktop

If nothing happens, download GitHub Desktop and try again.

Launching Xcode

If nothing happens, download Xcode and try again.

Launching Visual Studio Code

Your codespace will open once ready.

There was a problem preparing your codespace, please try again.

Latest commit

Git stats

Files

Failed to load latest commit information.

README.md

Решение краевой задачи при помощи языка Python.

Этот проект должны выполнить студенты 3-го курса ПМ Добровольский Георгий и Кривошеева Марина.

Кривошеева Марина — разработка вычислительной части (решение задачи методом конечных элементов). Добровольский Георгий — интерфейс программы, построение графика найденного решения.

Графический интерфейс реализован при помощи PyQt5 и Matplotlib . Соответственно, для запуска программы понадобятся эти библиотеки. Кроме того, в решении самой задачи используется библиотека Numpy .

Источник

Метод конечных разностей

Необходимо решить задачу с помощью метода конечных разностей на Python:

Моделирование квантовых систем: Метод конечных разностей

Реализовать процедуру solveFDM(mesh, potential, num_levels), которая находит уровни энергии и волновые функции в заданном потенциале методом конечных разностей. На входе: mesh = [x1, x2, . , xN] – сетка, заданная в виде массива точек; potential = [V1, V2, . , VN] – потенциал, заданный в точках сетки; num_levels – число состояний, которое нужно найти. Считать, что на краях отрезка (т.е. в точках x = 0 и x = L) волновая функция обращается в ноль.

Написать процедуры, которые строят график потенциала и волновых функций.

Протестировать процедуру на каком-нибудь потенциале (на ваш выбор).

f’’(xi) + V(xi)*f(xi) = E*f(xi) — это уравнение нужно решить с помощью метода конечных разностей, где f(xi) — волновая функция

Метод конечных элементов
Здраствуйте, мне нужно написать программную реализацию системы дифференциальных уравнений типа.

Метод конечных разностей
Задано дифференциальное уравнение у»+(3х+4)*у’+(х-5)*у=(х+2)^2 второго порядка с двухточечными.

Метод конечных разностей
Интерполяция методом конечных разностей на C++. Подскажите алгоритм. И необходимо сделатьввод.

Метод конечных разностей
Здравствуйте! Нужна помощь с методом конечных разностей. На этом сайте.

Так «помочь решить задачу» или «полностью решить задачу вместо вас»?
Если первое — задавайте вопросы, которые неясны. Желательно по одному.
Если второе — то вам в раздел «предлагаю работу».

Источник

Функции конечных разностей Python?

Я’искал в Numpy/Scipy модули, содержащие функции конечных разностей. Однако, самое близкое, что я нашел — это numpy.gradient() , который хорош для конечных разностей 1-го порядка точности 2-го порядка, но не очень, если вам нужны производные более высокого порядка или более точные методы. Я даже не нашел очень много специальных модулей для такого рода вещей; большинство людей, похоже, делают «накрутку» по мере необходимости. Поэтому мой вопрос заключается в том, знает ли кто-нибудь модули (либо часть Numpy/Scipy, либо сторонний модуль), которые специально посвящены методам конечных разностей высшего порядка (как по точности, так и по производной). У меня есть свой собственный код, над которым я работаю, но он в настоящее время довольно медленный, и я не собираюсь пытаться оптимизировать его, если уже есть что-то доступное.

Обратите внимание, что я говорю о конечных разностях, а не о производных. Я видел и scipy.misc.derivative() и Numdifftools, которые берут производную аналитической функции, которой у меня нет.

Один из способов сделать это быстро — свертка с производной гауссова ядра. Простым случаем является свертка вашего массива с [-1, 1] , что дает простую формулу конечных разностей. Кроме того, (f*g)’= f’*g = f*g’ , где * — свертка, так что в итоге вы получаете производную, свернутую с обычным гауссом, что, конечно, немного сгладит ваши данные, что можно минимизировать, выбрав наименьшее разумное ядро.

import numpy as np from scipy import ndimage import matplotlib.pyplot as plt #Data: x = np.linspace(0,2*np.pi,100) f = np.sin(x) + .02*(np.random.rand(100)-.5) #Normalization: dx = x[1] - x[0] # use np.diff(x) if x is not uniform dxdx = dx**2 #First derivatives: df = np.diff(f) / dx cf = np.convolve(f, [1,-1]) / dx gf = ndimage.gaussian_filter1d(f, sigma=1, order=1, mode='wrap') / dx #Second derivatives: ddf = np.diff(f, 2) / dxdx ccf = np.convolve(f, [1, -2, 1]) / dxdx ggf = ndimage.gaussian_filter1d(f, sigma=1, order=2, mode='wrap') / dxdx #Plotting: plt.figure() plt.plot(x, f, 'k', lw=2, label='original') plt.plot(x[:-1], df, 'r.', label='np.diff, 1') plt.plot(x, cf[:-1], 'r--', label='np.convolve, [1,-1]') plt.plot(x, gf, 'r', label='gaussian, 1') plt.plot(x[:-2], ddf, 'g.', label='np.diff, 2') plt.plot(x, ccf[:-2], 'g--', label='np.convolve, [1,-2,1]') plt.plot(x, ggf, 'g', label='gaussian, 2')

Поскольку вы упомянули np.gradient , я предположил, что у вас есть по крайней мере 2d массивы, поэтому следующее относится к ним: Это встроено в пакет scipy.ndimage , если вы хотите сделать это для ndarrays. Будьте осторожны, потому что, конечно, это не даст вам полный градиент, но я полагаю, что это произведение всех направлений. Надеюсь, кто-нибудь с лучшими знаниями выскажется.

from scipy import ndimage x = np.linspace(0,2*np.pi,100) sine = np.sin(x) im = sine * sine[. None] d1 = ndimage.gaussian_filter(im, sigma=5, order=1, mode='wrap') d2 = ndimage.gaussian_filter(im, sigma=5, order=2, mode='wrap') plt.figure() plt.subplot(131) plt.imshow(im) plt.title('original') plt.subplot(132) plt.imshow(d1) plt.title('first derivative') plt.subplot(133) plt.imshow(d2) plt.title('second derivative')

Использование gaussian_filter1d позволяет взять направленную производную вдоль определенной оси:

imx = im * x d2_0 = ndimage.gaussian_filter1d(imx, axis=0, sigma=5, order=2, mode='wrap') d2_1 = ndimage.gaussian_filter1d(imx, axis=1, sigma=5, order=2, mode='wrap') plt.figure() plt.subplot(131) plt.imshow(imx) plt.title('original') plt.subplot(132) plt.imshow(d2_0) plt.title('derivative along axis 0') plt.subplot(133) plt.imshow(d2_1) plt.title('along axis 1')

Первый набор результатов, приведенных выше, меня немного смущает (пики в оригинале отображаются как пики во второй производной, когда кривизна должна быть направлена вниз). Не изучив подробнее, как работает 2-я версия, я могу порекомендовать только 1-ю версию. Если вам нужна величина, просто сделайте что-то вроде:

d2_mag = np.sqrt(d2_0**2 + d2_1**2)

Источник