Решить дифференциальное уравнение методом Рунге-Кутта 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]);

@Moro1 · 17.01.2014, 21:23 **[ТС]**

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

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

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

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

Vy pytaetes obyasnit, eto yge zamechatelno!!!
T.e pravilno li ya ponimayu dq(1) - eto (prostymi slovami) konteiner v kotoryi popadyt znacheniya v kagdyi moment vremeni, posle resheniya yravneniya s pravoi storony (v moem slychae V-R*q(1)

? Esli eto tak, a q(1) - pervaya proizvodnaya, q(2)- vtoraya proizvodnaya i t.d., togda ne ponyatno, chto stroit matlab. Xotyaby vzyat dq(3)=V-R*q(2); -= reshenie dolgno byt nylevoe (ibo vtoraya proizvodnaya ot vremeni =0), a on mne vydaet kakyuty zagagyliny. neponyatno.

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
function tempq
[T,Q] = ode45(@loreq, [0,3], [0 2 1]);
plot(T,Q)
end
 
function dq = loreq(t,q) 
R = 2;
V=4;
dq=zeros(3,1)
dq(1)=V-R; 
dq(2)=V/R*q(1);
dq(3)=V-R*q(2);
end

@tanyabo · 17.01.2014, 23:57

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

kak programma opredelyaet to chto nygno

это решение задачи Коши, в которой говорится, что мы первую производную dy(1) обозначаем за y(2), вторую производную за y(3), третью - за y(4) и так до тех пор пока не достигнем порядка уравнения. Т е. например, если у тебя уравнение третьего порядка, то (это правило такое решения диффур, просто прими это

)так надо, запомни раз и все)
dy(1) = y(2)
dy(2) = y(3)
а dy(3) = .....само уравнение, в котором ты заменяешь, например, y' на y(2) а y'' на y(3)
ты же понимаешь, что dy(1) это y' ??
вроде так

я так решала все диффуры, может я не права, но как объяснить еще не знаю

@Moro1 · 18.01.2014, 00:46 **[ТС]**

tanyabo, Огромное спасибо, что откликнулись!!!
Однако... Несмотря на то, что уже так делал. Сделаю еще раз.
Вы хотите сказать, что для уравнения (закон Ома): V-R*d2q/dt2=0 получится кривая отмеченная красной стрелочкой? А для уравнения V-R*dq/dt=0, кривая отмеченная оранжевой стрелочкой? Или я недопонимаю чего-то?

@tanyabo · 18.01.2014, 02:38

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

я недопонимаю чего-то

или я не понимаю чего-то)) откуда у тебя 4 начальных условия и 4 производные в программе, если у тебя уравнение второго порядка

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

d2q/dt2

а чтобы понять что есть что на графике, построй зависимости не просто Q от Т, а зависимости каждой производной от T отдельно

@Moro1 · 18.01.2014, 03:15 **[ТС]**

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

или я не понимаю чего-то)) откуда у тебя 4 начальных условия и 4 производные в программе, если у тебя уравнение второго порядка
а чтобы понять что есть что на графике, построй зависимости не просто Q от Т, а зависимости каждой производной от T отдельно

Зависимость каждой производной от времени выполняется как: [T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options) ???
Или так:

Matlab M
1
2
3
4
function tempq
[T,Q] = ode45(@loreq, [0,3], [2 0 0 2]);
plot(T,Q(:,1),T,Q(:,2),T,Q(:,3),T,Q(:,4))
end

???

Зосима · 18.01.2014, 13:10

Moro1, ту систему из 4х уравнений можно трактовать, как систему двух функций a(t) и b(t):
a''(t) = V-R*a'(t); b''(t) = V-R*b'(t);
лишь из-за того, что начальные условия разные получаем разные графики

если ты решаешь закон Ома: d2q/dt2 = V/R; то у тебя получается система из двух уравнений и в результате ты поимеешь два графика: саму функцию q(t) и первую прозводную q'(t).

*подписывать кривые лучше всего легендой (на примере моей безумной системы 5го порядка)

Matlab M
1
2
plot(T,Y)
legend('y(t)','первая прозводная dy','d2y','d3y','d4y');

@Moro1 · 18.01.2014, 14:24 **[ТС]**

Все вроде разобрался. Всем спасибо.
Если работать по методу предложенному tanyabo, - то все графики неверны, не может так заряд себя вести (первая картинка).
Приложил картинку (№2) в которой все получается как надо!

@tanyabo · 18.01.2014, 14:34

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

все графики неверны

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

Зосима · 18.01.2014, 15:26

Moro1, tanyabo, та система с синусообразными графиками - это система из двух различных функций второго порядка!
равно как и вторая программа- это система из 3х функций первого порядка!!!

напиши просто формулу (а не прогрвмму) а я подскажу как это должно быть

@Moro1 · 18.01.2014, 16:48 **[ТС]**

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

Moro1, tanyabo, та система с синусообразными графиками - это система из двух различных функций второго порядка!
равно как и вторая программа- это система из 3х функций первого порядка!!!

напиши просто формулу (а не прогрвмму) а я подскажу как это должно быть

Так я не понял, я закодировал неправильно?
В общем, то что я пытаюсь сделать это не задача, как таковая. Я просто пытаюсь понять как работае Matlab.
Я хочу построить 3 графика (напримере закона Ома)!
Первый, зависимость q(t), второй, зависимость dq/dt(t), третий, зависимость I(V)

Зосима · 18.01.2014, 17:41

да, с системой напортачил чуток

нарисовать q(t) и dq/dt что есть I(t) - проблем нет

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
function tempq 
[T,Q] = ode45(@loreq, [0,3], [0 2]); 
plot(T,Q) 
legend('q(t)', 'I(t)' )
end 
 
function dq = loreq(t,q) 
R = 2; 
V=4; 
dq=zeros(2,1) 
dq(1)=q(2);
dq(2)=V/R;
end

а вот I(V) трудно изобразить, т.к. V=const

@Moro1 · 18.01.2014, 18:21 **[ТС]**

А вот так получается, что система строит 'q(t)', 'dqdt(t)', 'd2qdt(t)', 'd3qdt(t)'. Верно?

Matlab M
1
2
3
4
dq(1)=q(2);
dq(2)=q(3)
dq(3)=q(4)
dq(4)=V/R;

Зосима · 18.01.2014, 18:47

да, верно

это равносильно ОДУ: d4q/dt4 = V/R

@Moro1 · 18.01.2014, 19:57 **[ТС]**

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

а вот I(V) трудно изобразить, т.к. V=const

Правильно ли я понимаю, что если сделать V сделать плавающей величиной, то код придется переписывать полностью. Т.е конструкция:

Matlab M
1
2
3
[T,Q] = ode45(@loreq, [0,3], [0 2]); 
...................................................................
function dq = loreq(t,q)

уже не подойдет?

Зосима · 18.01.2014, 21:08

не, там чуток подправить

гляди:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function tempq 
[ T , Q ] = ode45 ( @loreq , [0, 3] , [0, 2] ) ; 
plot( T , Q, V(T), Q(:,2) ) 
legend( ' q(t) ' , ' I(t) ' , 'I(V)') 
end
 
function dq = loreq( t , q ) 
R = 2 ; 
dq = zeros( 2 , 1 ) ;
dq( 1 ) = q( 2 ) ; 
dq( 2 ) = V(t) / R ; 
end
 
function y = V(t) 
y = sin(2*pi*4*t).*(t<=1);
end

видишь, я только добавил функцию напряжения от времени

*только график покажи - любопытно

@Moro1 · 18.01.2014, 22:06 **[ТС]**

Немного подкорректировал программу, ибо при смеси q(t), dqdt и I(V) получается каша

) А при таком варианте: милое колесо иезекииля.
Вы не могли бы прокомментировать, что значит запись

Matlab M
1
(t<=1);

Зосима · 20.01.2014, 12:09

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

что значит запись
(t<=1);

Если ты построишь график V(t), то заметишь, что это синус, который после t=1 становится нулем, т.е. фактически радиоимпульс длительностью 1

Название: 0f3.png
Просмотров: 64

Размер: 1.5 Кб

А кружок получился потому как токи напряжение имеют сдвиг фазы pi/2: U(t)~sin(t), I(t)~cos(t), и получаем фигуру Лиссажу

вот как они будут рядом выглядеть:

Кликните здесь для просмотра всего текста

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
function temp
tk = 0:0.01:3;
[ T , Q ] = ode45 ( @loreq , tk , [0, 2] ) ; 
plotyy( T , Q(:,2), T, V(T)) 
legend(' I(t) ' , 'V(t)') 
end
 
function dq = loreq( t , q ) 
R = 2 ; 
dq = zeros( 2 , 1 ) ;
dq( 1 ) = q( 2 ) ; 
dq( 2 ) = V(t) / R ; 
end
 
function y = V(t) 
y = sin(2*pi*4*t).*(t<=1);
end

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

@Moro1 · 21.01.2014, 11:55 **[ТС]**

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

да, с системой напортачил чуток

нарисовать q(t) и dq/dt что есть I(t) - проблем нет

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
function tempq 
[T,Q] = ode45(@loreq, [0,3], [0 0]); 
plot(T,Q) 
legend('q(t)', 'I(t)' )
end 
 
function dq = loreq(t,q) 
R = 2; 
V=4; 
dq=zeros(2,1) 
dq(1)=q(2);
dq(2)=V/R;
end

а вот I(V) трудно изобразить, т.к. V=const

Великодушно меня извините. Я все с Matlabом бьюсь. Решил это уравнение вручную, приложил картинку. Все идентично, в программе и на бумаге. Почему, тогда разные результаты? Никак не пойму принцип расчета. софта.

@Moro1 · 21.01.2014, 18:19 **[ТС]**

Ya popytalsya reshit cherez funkciyu diff, no pochemyto nichego ne vyxodit.

Matlab M
1
2
3
4
5
6
7
8
syms x;
V=4;
R=2;
f=V/R;
f2=diff(f,2);
x=0:0.1:3;
plot(x,f2);
grid;

Добавлено через 20 минут
Vo tak eshe poproboval, no grafik pochemyto vyvodit kak q(t), xotya ykazal kak q''(t)

Matlab M
1
2
3
4
5
6
V=25;
R=2;
x=V/R;
t=1:1:20;
f = @(x) diff(f(x),2);
plot(t,f);

Зосима · 22.01.2014, 11:40

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

Решил это уравнение вручную, приложил картинку. Все идентично, в программе и на бумаге.

эх, где же идентично то?

систему не правильно написал: в программе d2q/dt2 = V/R, а на бумаге dq/dt = V/R

Добавлено через 8 минут
вот как его рисовать:

Matlab M
1
2
3
4
5
6
7
syms x;
V=4;
R=2;
Q = dsolve('Dy=V/R', 'y(0)=0');
Q = subs(Q)
ezplot(char(Q),[0 3]);
grid;

Новые блоги и статьи Все статьи Все блоги /
Запустил конкурс "тем и промптов для текстовых квестов созданных почти чисто ИИ" Adler 06.10.2026 Всем привет! За последние три-четыре дня я создал более 16 текстовых квестовых игр используя преимущественно по одному запросу к ИИ на игру. Мне так понравилось смотреть все ветки/ сцены во всех. . .	ИИ не может найти нужный язык в списке Supersumestria 05.10.2026 Я ему даю вот такое изображение и прошу найти и подчеркнуть немецкий язык. Возвращает он вот это: https:/ / i. **********/ vqBWLe2. png Нужную строчку в 3й колонке просто выдумал. . Это. . .	Новая последняя моя музыка в SUNO zorxor 05.10.2026 Здравствуйте, дорогие мои друзья! С большой радостью я хотел бы представить вам свою новую последнею музыку, которую сгенерировала мне по моей просьбе нейросеть SUNO. С уважением, zorxor. Это. . .	Nekobox - outbounds[0].transport: unknown transport type: raw damix 01.10.2026 Фикс ошибки Правым кликом по серверу -> отладочная информация -> edit Заменить "net": "raw", на "net": "tcp", Нажать кнопку reload.
Программный домашний кинотеатр russiannick 27.09.2026 Сподобился на программный домашний кинотеатр. В качестве ЯВУ по традиции выбрал js. В помощники взял Яндекс-Алису. Было создано три зала на разные интересы. исторические и ретро сериал Хичкок. . .	Беседа с ИИ о программистах, недопускающих к созданию и правке кода генеративные ИИ и причины этого zorxor 21.09.2026 Раньше я радовался или получал некоторые эмоции, пусть небольшие, но всё же, от самого процесса написания кода, рекомпиляции и запуска, видя постепенное развитие программы и прочее. А теперь лень. . .	Мобильное приложение ColorStep pavlinmavlin 17.09.2026 Реализовал приложение Красный, Зеленый, Синий в Unity3d + c#. Название изменил на ColorStep. Приложение прошло модерацию и теперь доступно для скачивания. Делал его сам, шаг за шагом — и вот,. . .	Запрет дублирования строк в табличной части Maks 13.09.2026 Реализация из решения ниже выполнена на нетиповом справочнике "Нормы ТО" с табличной часть "Виды ТО", разработанного в КА2, со следующими реквизитами: - ВидТО (СправочникСсылка. ВидыТО); - ВидГСМ. . .

@Moro1 0 / 0 / 0 Регистрация: 31.12.2013 Сообщений: 48
	18.01.2014, 00:46 [ТС]
	tanyabo, Огромное спасибо, что откликнулись!!! Однако... Несмотря на то, что уже так делал. Сделаю еще раз. Вы хотите сказать, что для уравнения (закон Ома): V-Rd2q/dt2=0 получится кривая отмеченная красной стрелочкой? А для уравнения V-Rdq/dt=0, кривая отмеченная оранжевой стрелочкой? Или я недопонимаю чего-то? Миниатюры 0

@Moro1 0 / 0 / 0 Регистрация: 31.12.2013 Сообщений: 48
	18.01.2014, 14:24 [ТС]
	Все вроде разобрался. Всем спасибо. Если работать по методу предложенному tanyabo, - то все графики неверны, не может так заряд себя вести (первая картинка). Приложил картинку (№2) в которой все получается как надо! Миниатюры 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	18.01.2014, 15:26
	Moro1, tanyabo, та система с синусообразными графиками - это система из двух различных функций второго порядка! равно как и вторая программа- это система из 3х функций первого порядка!!! напиши просто формулу (а не прогрвмму) а я подскажу как это должно быть 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	18.01.2014, 18:47
	да, верно это равносильно ОДУ: d4q/dt4 = V/R 1