Решение системы дифференциальных уравнений методом Эйлера и Рунге-Кутты 4 порядка - F

Алгоритмы: от сортировки пузырьком до численных методов

Модераторы: C_O_D_E, DeeJayC

Ответить
Turb_man
Сообщения: 9
Зарегистрирован: 30 июн 2017, 07:15

30 июн 2017, 07:23

Здравствуйте уважаемые форумчани.

Есть следующая задача: решить методом Эйлера и Рунге-Кутты 4 порядка следующие уравнение

dv/dt = - (sx + v)/tau, где x=x(t), v=v(t) - координата и скорость движения частицы, t- время движения, s, tau - произвольные параметры.

Начальные условия x(0)=x0, v(0)=0.
Так как x = dv/dt, то получается уравнение x"+x/tau+s/tau*x=0. Точное решение этого уравнения x = x0 * exp(-t/(2*tau))*cos(sqrt(s/tau-1/(4*tau^2)*t)).

Вроде бы я точное решение нашёл верно.

Теперь когда я решаю данное уравнение с помощью метода Эйлера и Рунге-Кутта 4 порядка, то для достаточно маленького шага численное решение довольно-таки значительно отличается, от точного, хотя при таком шаге и числе точек все численные и точные решения должны совпадать.

Вот код

program interpolation
integer i, j
real*8 x0, tau, x(1000), s, dt
character(6) file_name
c Запись точного решения в файл
file_name = 'myfile'
c Задание начальных условий
x0 = 5
tau = 3
s = 5
dt = 0.02
open(unit = 1, file = file_name)
do i = 1, 1000
c Аналитичсекое решение уравнения
x(i) = x0*exp(-(i-1)*dt/(2*tau))*cos(sqrt(s/tau-1/
&(4*tau**2)) * (i-1)*dt)
write(1,*) (i-1)*dt, x(i)
enddo
close(1)
c Вызов функций интерполяции
call rungi_kutt
call eler
end
c Метод Рунге-Кутты
subroutine rungi_kutt
c Использование внешней функции
external F
real*8 k1(2), k2(2), k3(2), k4(2), F
real*8 y(1000,2), dt
integer i, j
character(10) file_name
c Имя файла для записи данных
file_name = 'rk_file'
c Шаг по времени
dt = 0.02
y = 0.0
c Открытие файла для записи данных
open(unit = 2, file = file_name)
c Алгоритм Рунге-Кутты 4 порядка
do i = 1, 999, 1
y(1,1) = 5.0
y(1,2) = 0.0
k1(1) = y(i,2)*dt
k1(2) = F(y(i,1), y(i,2))*dt
k2(1) = (y(i,2) + k1(2)/2)*dt
k2(2) = F(y(i,1)+k1(1)/2, y(i,2)+k1(2)/2)*dt
k3(1) = (y(i,2) + k2(2)/2)*dt
k3(2) = F(y(i,1)+ k2(1)/2, y(i,2)+ k2(2)/2)*dt
k4(1) = (y(i,2) + k3(2))*dt
k4(2) = F(y(i,1)+1*k3(1),y(i,2)+k3(2))*dt
do j = 1,2
y(i + 1,j) = y(i,j) + (k1(j) + 2*k2(j) + 2*k3(j) + k4(j))/6
enddo
enddo
c Запись данных в файл rk_file
do i = 1, 1000
write(2,*) (i-1)*dt, y(i+5,1)
enddo
close(2)
end
c-----------------------------------------------------------------------
function F(x,v)
real*8 x,v,f,s,tau
tau = 3
s = 5
f = 0
f = (-s*x-v)/tau
return
end
c-----------------------------------------------------------------------
c Эйлера
subroutine eler
c Использование внешней функции
external F
real*8 tau, s, dt, F
real*8 z(1000,2)
integer i, j
character(10) file_name
c Имя файла для записи данных
file_name = 'eler_file'
dt = 0.02
z(1,1) = 5.0
z(1,2) = 0.0
c Открытие файла для записи данных
open(unit = 1, file = file_name)
c Алгоритм Эйлера
do i = 1, 999
z(i+1,1) = z(i,1) + dt * z(i,2)
z(i+1,2) = z(i,2) + F(z(i,1), z(i,2)) * dt
enddo
c Запись данных в файл eler_file
do i = 1, 1000
write(1,*) (i-1)*dt, z(i,1)
enddo
close(1)
end
c-------------------------------------------------------------------------------------------------------------
Я в численных методах не очень силён, но просмотрел все сообщения данного форума по этой теме, но не удалось найти решение. Также прошу не кидать ссылки на статьи из викепедии, я уже там был и всё перечитал.

Пожалуйста, если найдёте ошибку, то сообщите, где она. Сами алгоритмы интерполяции я проверял, вроде правильно.
Слива
Сообщения: 133
Зарегистрирован: 19 мар 2016, 10:15

30 июн 2017, 11:09

Это какой язык-то вообще?
Слива
Сообщения: 133
Зарегистрирован: 19 мар 2016, 10:15

30 июн 2017, 11:18

Цитата:"Так как x = dv/dt," - почему у Вас координата x это производная от скорости по времени. Вообще-то наоборот. Скорость - это производная от координаты x по времени. А ускорение - это вторая производная от координаты x по времени. Понятно?
Аватара пользователя
AiK
Сообщения: 2273
Зарегистрирован: 13 фев 2004, 18:14
Откуда: СПб
Контактная информация:

30 июн 2017, 12:34

Turb_man писал(а):Здравствуйте уважаемые форумчани.

Есть следующая задача: решить методом Эйлера и Рунге-Кутты 4 порядка следующие уравнение

dv/dt = - (sx + v)/tau, где x=x(t), v=v(t) - координата и скорость движения частицы, t- время движения, s, tau - произвольные параметры.

Начальные условия x(0)=x0, v(0)=0.
Так как x = dv/dt, то получается уравнение x"+x/tau+s/tau*x=0. Точное решение этого уравнения x = x0 * exp(-t/(2*tau))*cos(sqrt(s/tau-1/(4*tau^2)*t)).

Вроде бы я точное решение нашёл верно.
Это однородное дифференциальное уравнения второго порядка с постоянными коэффициентами. Точных решений у этого уравнения не 1, а 3, в зависимости от знака дискриминанта характеристического уравнения. Выбран только третий вариант и в нём явно не хватает упоминания синуса с тем же аргументом, что и косинус.
Даже самый дурацкий замысел можно воплотить мастерски
Turb_man
Сообщения: 9
Зарегистрирован: 30 июн 2017, 07:15

03 июл 2017, 06:41

Спасибо большое всем, кто откликнулся. Фазу я действительно забыл. Мне уже подсказали. Методы сами в целом правильные.
Turb_man
Сообщения: 9
Зарегистрирован: 30 июн 2017, 07:15

03 июл 2017, 06:45

Слива писал(а):Это какой язык-то вообще?
Это Fortran. Так как x = dv/dt - извините, здесь я просто описался.
Ответить