Решить дифференциальное уравнение методом Рунге-Кутта 4 порядка (ode45)

@Moro1 · Регистрация: 31.12.2013

Студворк — интернет-сервис помощи студентам

Добрый день!
Необходимо решить дифференциальное уравнение методом Рунге-Кутта 4 порядка.
Уравнение подготовил, а вот как реализовать его через функцию ode45, с построением графика не знаю.
Подскажите пожалуйста, какой код необходимо дописать?
Код программы:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
e=15;
V1=10;
t=(1:0.1:30);
V2=1;
f=10000;
omega=2*pi*f;
k=1.36*10-23;
Lk=0.000001;
R=1600;
T=222;
B=222;
Vc=0.92837;
R=1;
Vn=sqrt(4*k*R*T*B);
function q(2)=solverDE(t,y);
f=@(t,q)([q(2); V1+V2*sin*(omega*t)-Lk*q(2)-R*q(1)-Vn-Vc*sin((2*pi*q(1))/2*e)/Lk]);

@Галина Борисовн · 01.01.2014, 12:26

Я немного откорректировала вашу программу: добавила начальные условия; взяла модуль под корнем. Теперь работает.

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
clear; clc
e=15;
V1=10;
t=[1:0.01:3];
V2=1;
f=10000;
omega=2*pi*f;
k=1.36*10-23;
Lk=0.000001;
R=1600;
T=222;
B=222;
Vc=0.92837;
R=1;
Vn=sqrt(abs(4*k*R*T*B));
x0=[0 0.5];%Начальные условия
f=@(t,q) [q(2); V1+V2*sin(omega*t)-Lk*q(2)-R*q(1)-Vn-Vc*sin((2*pi*q(1))/2*e)/Lk]; 
[T,X]=ode45(f,t,x0);
plot(T,X(:,1))

@Moro1 · 15.01.2014, 14:47 **[ТС]**

original message

Izvinite chto translitom.
Pomogite razobratsya s ewe odnoi vesch'yu. Chto v matlabe est q(2) i chto est q(1) (15 stroka)? Pochemy sistema provodit rasschet ne imeya dannyx o etix peremennyx, prichem samoe interesnoe, sistema stroit iskluchitelno pri q1 i q2. q3 i t.d. - vydaet owibky. Xotel razobratsya podrobnee. Zaranee spasibo!

Помогите разобраться с еще одной вещью. Что в матлабе есть q(2) и что есть q(1) (15 строка)? Почему система проводит расчет не имея данных о этих переменных, причем самое интересное, система строит исключительно при q1 и q2. q3 и тд - выдает ошибку. Хотел разобраться подробнее. Заранее спасибо!

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
clear; clc
e=1.6*10-19;
V1=1000;
t=0:1*10-9:1*10-6;
f=1*10+9;
omega=2*pi*f;
k=1.36*10-23;
Lk=0.000001;
R1=1.3*10+6;
T=0.03;
B=1*10+9;
R=1;
V2=sqrt(abs(4*k*R*T*B));
x0=[0.10 20];
f=@(t,q) [q(2);(q(2))]; 
[T,X]=ode45(f,[0 1],x0);
plot(T,X(:,1))
grid on;

@Галина Борисовн · 15.01.2014, 14:53

f=@(t,q) [q(2);(q(2))];
Эта запись мне не понятна.

@Moro1 · 15.01.2014, 15:08 **[ТС]**

Сообщение от Галина Борисовн

f=@(t,q) [q(2);(q(2))];
Эта запись мне не понятна.

original message

Vot i ya pro to govoryu. Vse delo v tom, chto iznachalno mne nygno bylo zaprogrammirovat formuly (prilogil kartinky). Chto sobstvenno ya i sdelal. Odnako v processe voznikli trydnosti ya ne znal kak v matlabe zaprogrammirovat proizvodnyu (dq/dt i d^2/dt^2). Poprosil cheloveka podskazat, seichas vozmognosti svyazatsya s nim net. Na chto polychil otvet: q(1)=q, a q(2)=dq/dt. V tot moment somnenii ne vozniklo. Vrode rabotaet. A vot kogda stal proveryat polychennyi rezyltat doskonalno vyshla beliberda.

Вот и я про то говорю. Все дело в том, что изначально мне нужно было запрограммировать формулу (приложил картинку). Что собственно я и сделал. Однако в процессе возникли трудности я не знал как в матлабе запрограммировать производную (dq/dt и d^2/dt^2). Попросил человека подсказать, сейчас возможности связаться с ним нет. На что получил ответ: q(1)=q, a q(2)=dq/dt. В тот момент сомнений не возникло. Вроде работает. А вот когда стал проверять полученный результат досконально вышла билиберда.

@tanyabo · 15.01.2014, 16:09

Тебе непонятно почему q1=q а q2=dq/dt ?
Я решала методом рунге кутта, задай вопрос поконкретнее)
Тебе точно подходит ode45?

@Moro1 · 15.01.2014, 21:26 **[ТС]**

Сообщение от tanyabo

Тебе непонятно почему q1=q а q2=dq/dt ?
Я решала методом рунге кутта, задай вопрос поконкретнее)
Тебе точно подходит ode45?

На счет точности метода не знаю, почитал разные сайты, почему то решил что подойдет ode45. А что еще можно использовать? Мне непонятно откуда берется значения в матлабе для q(1) и q(2) если я иx не задавал.

Добавлено через 4 часа 52 минуты

Подскажите пожалуйста. Потиxоньку разбираюсь в матлабе. Пытаюсь решить диф. уравнение: V=R*dq/dt. Что делаю не так?

Matlab M
1
2
3
4
5
6
function z = loreq(t,q) 
global R; R = 10;
 z(1) = 560-R*q(1);
 
[T,Q] = ode45(@loreq, [0,3], [1])
plot(T,Q)

@bani-8 · 16.01.2014, 00:22

Сообщение от Moro1

podoidet ode45

решать можно с помощью функций описанных, например, здесь:
http://www.mathworks.com/help/... ode45.html

Когда я решала свою задачу, у меня не было сказано, что система жесткая, а на деле оказалось так, долго думала в чем проблема и почему ode45 не подходит, заменила на ode15s и все на свои места встало
А в твоем 7 посте у тебя диффура первого порядка
может стоит записать сперва в виде Коши?

Зосима · 16.01.2014, 16:50

Moro1, нужно расписать красивенько функцию системы

Matlab M
1
2
3
4
5
6
7
8
9
10
function dqdt = diffsys(t, q)
global e V1 V2 omega k Lk R T B Vc Vn
% обозначения:
% q(1) -> q(t) (искомая функция)
% q(2) -> dq/dt (первая производная)
% dqdt(1,1) -> dq/dt (первая производная)
% dqdt(2,1) -> d2q/dt2 (вторая производная)
dqdt(1,1) = q(2);
dqdt(2,1) = V1 + V2*sin(omega*t) - Lk*q(2) - R*q(1) - Vn - Vc*sin((2*pi*q(1))/2*e)/Lk; 
end

и программка рассчета:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
clear, clc
global e V1 V2 omega k Lk R T B Vc Vn
e=1.6*10^-19;  % это заряд электрона!!!
V1=10;
V2=1;
f=10000;
omega=2*pi*f;
k=1.36*10^-23; % здесь забыли степень!!!
Lk=0.000001;
R=1600;
T=222;
B=222;
Vc=0.92837;
Vn=sqrt(abs(4*k*R*T*B));
 
t = [0 3]; % диапазон времени
q0 = [0 0.5]; %Начальные условия
[T,X] = ode45('diffsys',t,q0);
plot(T,X(:,1))

В результате получаем убывающую чуть промодулированную синусоиду:

Решить дифференциальное уравнение методом Рунге-Кутта 4 порядка (ode45)

Сообщение от Moro1

Что я делаю не так?

наверное просто неверно составляешь файл-функцию

вначале должна идти основная ф-ция, а потом ф-ция системы:

Matlab M
1
2
3
4
5
6
7
8
9
function temp % основная ф-ция
[T,Q] = ode45(@loreq, [0,3], [1]);
plot(T,Q)
end
 
function z = loreq(t,q) % функция системы
R = 10;
z(1) = 560-R*q(1);
end

*это все один файл temp.m

Зосима · 16.01.2014, 17:42

Сообщение от Moro1

Мне не понятно откуда берутся значения в матлабе для q(1) и q(2) если я их не задавал.

Как же не задавал? Задавал!

ты задаешь начальные условия x0 и формулы для производный по которым шаг за шагом высчитываются новые значения функции, и низших производных!

В твоем случае диффура имеет второй порядок, т.к. фигурирует вторая производная. поэтому в качестве нач. условий нужно задать два значения: саму функцию q(1) (нулевая производная) и первую производную q(2) из которых уже будет рассчитана вторая производная, а по ней -значения ф-ции и первой производной в следующей точке!

и так по шагам двигаемся до конца временного интервала. Смекаешь?

Добавлено через 22 минуты

Сообщение от Moro1

Пытаюсь решить диф. уравнение: V=R*dq/dt

в этом случае система будет выглядеть так:
V = R*dq/dt;
dq/dt = V/R;

Matlab M
1
2
3
4
5
6
7
8
9
function temp % основная ф-ция
[T,Q] = ode45(@loreq, [0,3], [1]);
plot(T,Q)
end
 
function z = loreq(t,q) % функция системы
R = 10;
z(1) = 560/R;
end

и получаем прямую

@Moro1 · 16.01.2014, 18:14 **[ТС]**

Сообщение от Зосима

Как же не задавал? Задавал!

ты задаешь начальные условия x0 и формулы для производный по которым шаг за шагом высчитываются новые значения функции, и низших производных!

В твоем случае диффура имеет второй порядок, т.к. фигурирует вторая производная. поэтому в качестве нач. условий нужно задать два значения: саму функцию q(1) (нулевая производная) и первую производную q(2) из которых уже будет рассчитана вторая производная, а по ней -значения ф-ции и первой производной в следующей точке!

и так по шагам двигаемся до конца временного интервала. Смекаешь?

Добавлено через 22 минуты

в этом случае система будет выглядеть так:
V = R*dq/dt;
dq/dt = V/R;

Matlab M
1
2
3
4
5
6
7
8
9
function temp % основная ф-ция
[T,Q] = ode45(@loreq, [0,3], [1]);
plot(T,Q)
end
 
function z = loreq(t,q) % функция системы
R = 10;
z(1) = 560/R;
end

и получаем прямую

Mogno ytochnit.
pochemy z(1) = 560/R; esli, sydya po Vashim zapisyam vyshe
%q(1) -> q(t) (искомая функция)
%q(2) -> dq/dt (первая производная)
Po logike dolgno byt z(2) Chto to ya zapytalsya.

Зосима · 16.01.2014, 18:43

Нужно различать входные данные и выходные!

Входные - это t и q, а выходные тут z

на вход подаем t и q и получаем в ответ z

Добавлено через 8 минут
давай рассмотрим еще парочку наглядных диффур

Можешь записать систему для:

dy - exp(-t) = 5*y;

@Moro1 · 16.01.2014, 18:55 **[ТС]**

Сообщение от Зосима

Нужно различать входные данные и выходные!

Входные - это t и q, а выходные тут z

на вход подаем t и q и получаем в ответ z

Добавлено через 8 минут
давай рассмотрим еще парочку наглядных диффур

Можешь записать систему для:

dy - exp(-t) = 5*y;

Matlab M
1
2
3
4
5
6
7
8
9
function sad
[T,Q] = ode45(@loreq, [0,3], [1]);
plot(T,Q)
end
 
function dy = loreq(t,q) 
y = 10;
dy(1) = 5*y+exp(-t);
end

Зосима · 16.01.2014, 19:01

я подразумевал, что мы ищем не q(t), а y(t)

т.е. должно быть:

Matlab M
6
7
8
function dy = loreq(t,y) 
dy(1) = 5*y+exp(-t);
end

Но это не суть

теперь давай возьмем примерчик повеселее

d5y/dt5 + pi* d3y/dt3 - exp( d2y/dt2 - y ) + sin(t + y) = 0;

*ищем также y(t)

@Moro1 · 16.01.2014, 19:16 **[ТС]**

Я извиняюсь, пример слишком веселый

Что я растерялся...

Зосима · 16.01.2014, 19:30

тогда помогаю

Matlab M
1
2
3
4
5
6
7
function dy = loreq ( t , y )
dy = zeros(5,1); % создаем вектор столбец
dy(1)=y(2); 
dy(2)=y(3);
dy(3)=y(4);
dy(4)=y(5);
dy(5) = - pi* y(4) + exp( y(3) - y(1) ) - sin(t + y(1)) ;

@Moro1 · 16.01.2014, 20:48 **[ТС]**

Daite ewe kakoi-nibyd primer. Vrode ponyatno. A vrode... vot tak porobyi drygomy obyasnit - ne obyasnyu.

Добавлено через 48 минут
I ewe vopros. Chto est levaya edinica i chto est pravayu edinica (1,1; 2,1). Mne by kakyyunibyd ssylky na knigy ili material na saite, gde pro eto mogno prochitat. Spasibo zaranee.
% dqdt(1,1) -> dq/dt (первая производная)
% dqdt(2,1) -> d2q/dt2 (вторая производная)

I proverte pogalyista tak li ya zapisal: zelenaya liniya zavisimost dq/dt(t). Tak?

Matlab M
1
2
3
4
5
6
7
function z = loreq(t,q) 
R = 10;
V=30;
z=zeros(2,1);
z(1) = V/R;
z(2)=V-R*q(2);
end

Зосима · 17.01.2014, 16:00

Смотри

все матлабовские решатели odeXX требуют, чтобы функция системы возвращала вектор-столбец! Однако сделать столбец из набора данных (значений высших производных) можно разными стопобами

в хелпе для этого вначале создают столбец нулей нужной длинны:

Matlab M
dy = zeros(5,1); % 5 строк и 1 столбец

В этом же случае, я решил сократить код (ибо я очень ленивый

) и указать расположение элементов явно:
dqdt(1,1) - элемент в первой строке и в первом столбце
dqdt(2,1) - элемент во второй строке и в первом столбце
(это стандартная нумерация элементов в матлабе) поэтому в результате получаем массив dqdt в котором две строки и один столбец - т.е. именно то, что нужно!

Кроме того, сделать столбец можно еще и так, как это предлагала ГалинаБорисовна:

Matlab M
1
2
3
4
dqdt = [q(2)
      Vn + .... ];
% или
dqdt = [q(2);   Vn + .... ];

Тут создается массив (т.к. квадратные скобки), а элементы "укладываются" в столбец: так как впервом варианте идет перенос на новую строку, а во втором - точкозапятые - оба эти раздлителя говорят матлабу, что переходим на новую строку и получаем столбец. (см. тут)
Ясненько?

@Moro1 · 17.01.2014, 16:27 **[ТС]**

Prosti menya za moyu tverdolobost (mogete ne otvechat konechno, Vy i tak yge mnogo sdelali), kak togda programma opredelyaet stepeni proizvvodnye, chto est v skobkax: eto indeksy, ili ge eto stepeni? Ya prosto ne ponimayu kakaya svyaz mezdy y(3) i d3y/dt3, kak programma opredelyaet to chto nygno.
dy(1)=y(2);
dy(2)=y(3);
dy(3)=y(4);
dy(4)=y(5);
dy(5) = - pi* y(4) + exp( y(3) - y(1) ) - sin(t + y(1)) ;

Зосима · 17.01.2014, 16:33

тут dy(i) и y(i) - это ничто иное как просто элементы массива, индексы!

просто мы их сами обозначаем так, как нам удобно, т.е. номер производной равен или подобен номеру элемента

*просто я не преподаватель и не знаю, как доходчиво объяснить

Новые блоги и статьи Все статьи Все блоги /
Access VikBal 11.12.2025 Помогите пожалуйста !! Как объединить 2 одинаковые БД Access с разными данными.	Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro	Музыка, написанная Искусственным Интеллектом volvo 04.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .	От async/await к виртуальным потокам в Python IndentationError 23.11.2025 Армин Ронахер поставил под сомнение async/ await. Создатель Flask заявляет: цветные функции - провал, виртуальные потоки - решение. Не threading-динозавры, а новое поколение лёгких потоков. Откат?. . .	Поиск "дружественных имён" СОМ портов Argus19 22.11.2025 Поиск "дружественных имён" СОМ портов На странице: https:/ / norseev. ru/ 2018/ 01/ 04/ comportlist_windows/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .
Сколько Государство потратило денег на меня, обеспечивая инсулином. Programma_Boinc 20.11.2025 Сколько Государство потратило денег на меня, обеспечивая инсулином. Вот решила сделать интересный приблизительный подсчет, сколько государство потратило на меня денег на покупку инсулинов. . . .	Ломающие изменения в C#.NStar Alpha Etyuhibosecyu 20.11.2025 Уже можно не только тестировать, но и пользоваться C#. NStar - писать оконные приложения, содержащие надписи, кнопки, текстовые поля и даже изображения, например, моя игра "Три в ряд" написана на этом. . .	Мысли в слух kumehtar 18.11.2025 Кстати, совсем недавно имел разговор на тему медитаций с людьми. И обнаружил, что они вообще не понимают что такое медитация и зачем она нужна. Самые базовые вещи. Для них это - когда просто люди. . .	Создание Single Page Application на фреймах krapotkin 16.11.2025 Статья исключительно для начинающих. Подходы оригинальностью не блещут. В век Веб все очень привыкли к дизайну Single-Page-Application . Быстренько разберем подход "на фреймах". Мы делаем одну. . .	Фото: Daniel Greenwood kumehtar 13.11.2025

@Галина Борисовн 2835 / 2132 / 87 Регистрация: 02.05.2010 Сообщений: 3,194
	15.01.2014, 14:53
	f=@(t,q) [q(2);(q(2))]; Эта запись мне не понятна. 0

@tanyabo 3 / 3 / 0 Регистрация: 10.03.2013 Сообщений: 96
	15.01.2014, 16:09
	Тебе непонятно почему q1=q а q2=dq/dt ? Я решала методом рунге кутта, задай вопрос поконкретнее) Тебе точно подходит ode45? 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	16.01.2014, 18:43
	Нужно различать входные данные и выходные! Входные - это t и q, а выходные тут z на вход подаем t и q и получаем в ответ z Добавлено через 8 минут давай рассмотрим еще парочку наглядных диффур Можешь записать систему для: *dy - exp(-t) = 5y;** 0

@Moro1 0 / 0 / 0 Регистрация: 31.12.2013 Сообщений: 48
	16.01.2014, 19:16 [ТС]
	Я извиняюсь, пример слишком веселый Что я растерялся... 0

@Moro1 0 / 0 / 0 Регистрация: 31.12.2013 Сообщений: 48
	17.01.2014, 16:27 [ТС]
	Prosti menya za moyu tverdolobost (mogete ne otvechat konechno, Vy i tak yge mnogo sdelali), kak togda programma opredelyaet stepeni proizvvodnye, chto est v skobkax: eto indeksy, ili ge eto stepeni? Ya prosto ne ponimayu kakaya svyaz mezdy y(3) i d3y/dt3, kak programma opredelyaet to chto nygno. dy(1)=y(2); dy(2)=y(3); dy(3)=y(4); dy(4)=y(5); dy(5) = - pi* y(4) + exp( y(3) - y(1) ) - sin(t + y(1)) ; 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	17.01.2014, 16:33
	тут dy(i) и y(i) - это ничто иное как просто элементы массива, индексы! просто мы их сами обозначаем так, как нам удобно, т.е. номер производной равен или подобен номеру элемента *просто я не преподаватель и не знаю, как доходчиво объяснить 1

Решить дифференциальное уравнение методом Рунге-Кутта 4 порядка (ode45)

Решение