Решить дифференциальное уравнение методом Рунге-Кутта 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 · 05.02.2014, 17:51 **[ТС]**

Vopros. Pochemy v dannom slychae v funkcii I(V) polychaetsya vozrastayushaya pryamaya bez kakix-libo iskagenii, kotoraya daet 'sin'?
I vtoroe. Kak sdelat peremennyu 'q(1)' (na kartinke podrisoval dlya lychshego ponimaniya) zavisimyu ot 'function dq = zod(t,q)' v drygoi funkcii. Ili tak nelzya sdelat?

Matlab M
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
function tempqwertru
[T, Q] = ode45(@zod, [0,1e-12],[0 0]);
 
for i = 1:length(T)
    DQ(i,:) = zod(T(i),Q(i,:));
end
subplot(4,1,1); plot(T,Q(:,1));
legend('q(t)')
subplot(4,1,2); plot(T,DQ(:,1));
legend('dqdt(t)')
subplot(4,1,3); plot(T,DQ(:,2));
legend('d2qdt2(t)')
subplot(4,1,4); plot(T,V(T));
legend('I(V)')
end
 
function dq = zod(t,q)
L=1*10^-7; 
R=1*10^6;
dq=zeros(2,1);
dq(1)=q(2);
dq(2)=1/L*(V(t)-R*q(2));
end
 
function y = V(t) 
Va=1e-7;
Vc=1*10^-6;
y=(Va*t+Vc*sin((pi*1.6*10^-19)/1.6*10^-19));
end

@Moro1 · 05.02.2014, 19:12 **[ТС]**

Зосима, vy v proshlyi raz mne pokazali kak stroit zavisimost I(V(t)). Mne kagetsya, chto eto stroitsya kak to podrygomy, ibo plot(T,V(T)) - zavisimost vremne ot napryageniya po logike, ili ne tak? I odnovremenno ranshe yravneniya menee slognye reshala pravilno... V dobavok net literatyry kotoraya opisyvaet podobnye slychai...

Matlab M
1
2
3
4
5
6
7
8
subplot(4,1,4); plot(T,V(T));
legend('I(V)')
 
function y = V(t) 
Va=1e-7;
Vc=1*10^-6;
y=(Va*t+Vc*sin((pi*1.6*10^-19)/1.6*10^-19));
end

Добавлено через 10 минут
Moget I(V(t)) stroitsya kak plot(Q(:,2),V(T)); ?

Добавлено через 47 минут
Poka vopros odin i glavnyi. Kak sdelat peremennyu 'q(1)' zavisimyu ot 'function dq = zod(t,q)' v drygoi funkcii. Ili tak nelzya sdelat?

Matlab M
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
function tempqwertru
[T, Q] = ode45(@zod, [0,5e-13],[0 0]);
 
for i = 1:length(T)
    DQ(i,:) = zod(T(i),Q(i,:));
end
subplot(4,1,1); plot(T,Q(:,1))
subplot(4,1,2); plot(T,DQ(:,1))
subplot(4,1,3); plot(T,DQ(:,2))
subplot(4,1,4); plot(V(T),DQ(:,1))
legend('dqdt(V(t))')
end
 
 
function dq = zod(t,q)
L=1*10^-7; 
R=1*10^6;
dq=zeros(2,1);
dq(1)=q(2);
dq(2)=1/L*(V(t)-R*q(2));
end
 
function y=V(t);
Va=1*10^10;
Vc=1*10^-4;
y=(Va*t+Vc*sin(pi*q(1))/1.6*10^-19));
end

Зосима · 05.02.2014, 22:39

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

Moget I(V(t)) stroitsya kak plot(Q(:,2),V(T)); ?

Посмотри описание функции plot: http://www.exponenta.ru/soft/m... 0/plot.asp и все станет ясно

А вот это

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

Kak sdelat peremennyu 'q(1)' zavisimyu ot 'function dq = zod(t,q)' v drygoi funkcii

я попросту не понял

можешь простым математическим языком описать, что нужно?

*слух, а у тебя совсем не получается использовать экранную клавиатуру, бо транслит уже замучал

Добавлено через 2 минуты
Так, чтоли?

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
function dq = zod(t,q)
L=1*10^-7; 
R=1*10^6;
 
Va=1e-7;
Vc=1*10^-6;
V = Va*t+Vc*sin( pi*q(1)/1.6*10^-19 );
 
dq=zeros(2,1);
dq(1)=q(2);
dq(2)=1/L*(V-R*q(2));
end

@Moro1 · 05.02.2014, 23:08 **[ТС]**

С командами plot я давненько разобрался, просто почему то пользователям не дают удалять свои сообщения.
Теперь про математический язык, все что нужно на программном уровне есть под сполером.

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

Matlab M
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
function tempqwertru
[T, Q] = ode45(@zod, [0,5e-13],[0 0]);
 
for i = 1:length(T)
    DQ(i,:) = zod(T(i),Q(i,:));
end
subplot(4,1,1); plot(T,Q(:,1))
subplot(4,1,2); plot(T,DQ(:,1))
subplot(4,1,3); plot(T,DQ(:,2))
subplot(4,1,4); plot(V(T),DQ(:,1))
legend('dqdt(V(t))')
end
 
function dq = zod(t,q)
L=1*10^-7; 
R=1*10^6;
dq=zeros(2,1);
dq(1)=q(2);
dq(2)=1/L*(V(t)-R*q(2));
end
 
function y=V(t);
Va=1*10^10;
Vc=1*10^-4;
y=(Va*t+Vc*sin(pi*q(1))/1.6*10^-19));
end

Как сделать так, чтобы в 26 строке переменная q(1) бралась из функции выше, т.е с 14-20 строки. Ясно, что требуется?
О транслите: так претензий не было со стороны участников, надо было раньше сказать

)

Добавлено через 6 минут
При этом не меняя конструкцию кода.

Зосима · 05.02.2014, 23:11

Moro1, а как тебе тот вариант, что я предложил - считать V внутри zod?

Потом, когда ты получишь решение ОДУ, можно рассчитать и V(t), и V(I)

Не по теме:

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

О транслите: так претензий не было со стороны участников, надо было раньше сказать

виноват :-[ неважный с меня мойдодыр - буду на досуге твои творения переводить на литературный язык :pardon:

@Moro1 · 05.02.2014, 23:40 **[ТС]**

Зосима, Твой вариант, меня засмущал, как красну девицу

В прежней конструкции я разобрался, здесь опять темный лес, там было понятно "что за что в ответе". Тут же не ясно как выдернуть строчку 23 на ось абсцисс. С ординатой понятно...

Matlab M
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
function tempqwertru
[T, Q] = ode45(@zod, [0,1e-12],[0 0]);
 
for i = 1:length(T)
    DQ(i,:) = zod(T(i),Q(i,:));
end
subplot(4,1,1); plot(T,Q(:,1))
legend('dqdt(q(t))')
subplot(4,1,2); plot(T,DQ(:,1))
legend('dqdt(dqdt(t))')
subplot(4,1,3); plot(T,DQ(:,2))
legend('dqdt(d2qdt2(t))')
subplot(4,1,4); plot(V,DQ(:,1))
legend('dqdt(I(V))')
end
 
function dq = zod(t,q)
L=1*10^-7; 
R=1*10^6;
Va=1e-7;
Vc=1*10^-6;
 
V = (Va*t+Vc*sin(pi*q(1)/1.6*10^-19));
 
dq=zeros(2,1);
dq(1)=q(2);
dq(2)=1/L*(V-R*q(2));
end

Зосима · 06.02.2014, 00:14

ну тогда можно попробовать сделать так

Matlab M
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
function tempqwertru
[T, Q] = ode45(@zod, [0,5e-13],[0 0]);
 
for i = 1:length(T)
    DQ(i,:) = zod(T(i),Q(i,:));
    U(i) = V(T(i), Q(i,1));
end
subplot(4,1,1); plot(T,Q(:,1))
subplot(4,1,2); plot(T,DQ(:,1))
subplot(4,1,3); plot(T,DQ(:,2))
subplot(4,1,4); plot(U,DQ(:,1))
legend('dqdt(V(t))')
end
 
function dq = zod(t,q)
L=1*10^-7; 
R=1*10^6;
dq=zeros(2,1);
dq(1)=q(2);
dq(2)=1/L*(V(t,q(1))-R*q(2));
end
 
function y=V(t,q);
Va=1*10^10;
Vc=1*10^-4;
y=(Va*t+Vc*sin(pi*q)/1.6*10^-19));
end

@Moro1 · 06.02.2014, 00:26 **[ТС]**

И вопрос, почему матлаб строит зависимость I{V(t} как линейно возрастающую прямую, когда участвует эта штука sin((pi*q)/1.6*10^-19))?

Зосима · 06.02.2014, 11:03

Moro1, ну может это тебя натолкнет на размышления: при малых x => sin(x)~x

@Moro1 · 06.02.2014, 11:25 **[ТС]**

Зосима, Натолкнуло бы, при условии, что я не знал бы конечного результат. На выходе должно быть, что-то вроде ступенек, как раз таки за счет синуса и при таких исходных данных.
-----------------------

Matlab M
1
U(i) = V(T(i), Q(i,1));

- эта строка считает тоже производную? Как и эта:

Matlab M
1
DQ(i,:) = zod(T(i),Q(i,:));

?

Зосима · 06.02.2014, 11:36

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

эта строка считает тоже производную?

нет, она считает напряжение от времени и заряда

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

и еще если ток I~a*sin(w*t) и напряжение U~b*sin(w*t) то таки да - получим прямую

@Moro1 · 06.02.2014, 11:37 **[ТС]**

Matlab M
1
2
3
4
5
function y=V(t,q);
Va=1*10^10;
Vc=1*10^-4;
y=(Va*t+Vc*sin(pi*q)/1.6*10^-19));
end

Где 'y' есть dq/dt, верно?
--------
Нет, не верно

)

Зосима · 06.02.2014, 11:39

нет, конечно! y - это просто результат, который вернет функция/

@Moro1 · 17.02.2014, 19:17 **[ТС]**

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
function zaq
[T, Q] = ode45(@zod,[0,1e-15],[0 0]);
 
for i = 1:length(T)
    DQ(i,:)=zod(T(i),Q(i,:));
    U(i)=V(T(i),Q(i,1));
end
 plot(U,DQ(:,1))
legend('dqdt(V(t))')
end
 
function dq = zod(t,q)
L=1*10^-7; 
R=1*10^-6;
dq=zeros(2,1);
dq(1)=q(2);
dq(2)=1/L*(V(t,q(1))-R*q(2));
end
 
function y=V(t,q)
Va=1*10^-7;
Vc=1*10^-8;
y=(Va*t+Vc*sin(pi*q)/1.6*10^-19);
end

Скажите можно ли как то еще взят q из функции выше? Просто я никак не могу получит нужный мне результат... Почемуто даже изменение порядка в 22 строке совсем не меняет график, k primery s Vc=1*10^-8; k Vc=1*10^8; - результат одинаков

Зосима · 17.02.2014, 20:35

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

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

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

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

Если не сложно обьясните как программа понимает, что q берется из функции выше, а не отделная переменная.

Зосима · 17.02.2014, 20:56

ну как же!

в 17й строке мы обращаемся к функции V и отправляем в нее время t и заряд q(1)

@Moro1 · 18.02.2014, 20:56 **[ТС]**

Зосима, Понял в каком направлении нужно двигаться, и почему не получаеться нужныи результат, скажите как выполнить автоматическую очистку памяти?, команда clear all, почемуто не работает.
Картинки приложил, очевидно, что при отсутствии sin колебании не будет, a na grafikax net izmeneni sovsem.

Зосима · 20.02.2014, 11:00

А как ты понял, что clear не работает?

Синус возникает при решении ДУ второго порядка. Если напряжение постоянное, то зависимость I(U) тоже будет синусовое, если напряжение также имеет синусовость, то зависимость I(U) будет представлять фигуру Лиссажу.

@Moro1 · 20.02.2014, 13:05 **[ТС]**

Команда clear работает, но очищает все данные, после которого ничего не строится, выводится координатная ось и только.А хотелось бы чтобы программа очищала старые данные при каждом новом запуске.
Если напряжение постоянное, то зависимость I(U) тоже будет синусовое, если напряжение также имеет синусовость, то зависимость I(U) будет представлять фигуру Лиссажу.
Во вторых, по напряжению (на первой картинке), синуса - НЕТ!!! Обратите внимание, когда 0 умножается на синус (специально обвел линией), ну никак не может последний остаться. Но даже и это не главное, сравните два графика, они идентичны абсолютно. Разница должна же быть, Вы сами это подтвердили.

Новые блоги и статьи Все статьи Все блоги /
ИИ не может найти нужный язык в списке 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, со следующими реквизитами: - ВидТО (СправочникСсылка. ВидыТО); - ВидГСМ. . .	Скрипты Tampermonkey для CyberForum, ChatGPT, Claude и пр. Jin X 06.09.2026 Скрипты Tampermonkey для CyberForum, ChatGPT, Claude и пр. Работая с форумом и нейросетями в браузере часто хочется что-то подкорректировать или добавить какого-то функционала. Ниже прикреплён. . .

@Moro1 0 / 0 / 0 Регистрация: 31.12.2013 Сообщений: 48
	06.02.2014, 00:26 [ТС]
	И вопрос, почему матлаб строит зависимость I{V(t} как линейно возрастающую прямую, когда участвует эта штука *sin((piq)/1.610^-19))*? 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	06.02.2014, 11:03
	Moro1, ну может это тебя натолкнет на размышления: при малых x => sin(x)~x 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	06.02.2014, 11:39
	нет, конечно! y - это просто результат, который вернет функция/ 1

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	17.02.2014, 20:35
	дело в том, что логически программа составлена верно, я не знаю что тебе посоветовать 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	17.02.2014, 20:56
	ну как же! в 17й строке мы обращаемся к функции V и отправляем в нее время t и заряд q(1) 0

@Moro1 0 / 0 / 0 Регистрация: 31.12.2013 Сообщений: 48
	18.02.2014, 20:56 [ТС]
	Зосима, Понял в каком направлении нужно двигаться, и почему не получаеться нужныи результат, скажите как выполнить автоматическую очистку памяти?, команда clear all, почемуто не работает. Картинки приложил, очевидно, что при отсутствии sin колебании не будет, a na grafikax net izmeneni sovsem. Миниатюры 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	20.02.2014, 11:00
	А как ты понял, что clear не работает? Синус возникает при решении ДУ второго порядка. Если напряжение постоянное, то зависимость I(U) тоже будет синусовое, если напряжение также имеет синусовость, то зависимость I(U) будет представлять фигуру Лиссажу. 0

@Moro1 0 / 0 / 0 Регистрация: 31.12.2013 Сообщений: 48
	20.02.2014, 13:05 [ТС]
	Команда clear работает, но очищает все данные, после которого ничего не строится, выводится координатная ось и только.А хотелось бы чтобы программа очищала старые данные при каждом новом запуске. Если напряжение постоянное, то зависимость I(U) тоже будет синусовое, если напряжение также имеет синусовость, то зависимость I(U) будет представлять фигуру Лиссажу. Во вторых, по напряжению (на первой картинке), синуса - НЕТ!!! Обратите внимание, когда 0 умножается на синус (специально обвел линией), ну никак не может последний остаться. Но даже и это не главное, сравните два графика, они идентичны абсолютно. Разница должна же быть, Вы сами это подтвердили. 0