Решение системы дифф. уравнений и построение графиков

@Clarisse · Регистрация: 02.11.2012

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

Добрый вечер.Столкнулась с проблемой решения системы дифф. уравнений и построения графиков в Матлаб для своей курсовой работы. Системы из двух уравнений(две функции) решает на ура,теперь же нужно решить такую систему:

Code
1
2
3
4
5
6
7
8
a1*Q*v1 + a2*P*v2 + a3*P*v3 + a4*P*v4 - v1*(a1*v1+a2*v2+a3*v3+a4*v4)/(v1+v2+v3+v4)
a1*P*v1 + a2*Q*v2 + a6*P*v6 + a7*P*v7 - v2*(a1*v1+a2*v2+a6*v6+a7*v7)/(v1+v2+v6+v7)
a1*Q*v1 + a3*Q*v3 + a5*P*v5 + a6*P*v6 - v3*(a1*v1+a3*v3+a5*v5+a6*v6)/(v1+v3+v5+v6)
a1*Q*v1 + a4*Q*v4 + a5*P*v5 + a7*P*v7 - v4*(a1*v1+a4*v4+a5*v5+a7*v7)/(v1+v4+v5+v7)
a3*Q*v3 + a4*P*v4 + a5*Q*v5 + a8*P*v8 - v5*(a3*v3+a4*v4+a5*v5+a8*v8)/(v3+v4+v5+v8)
a2*Q*v2 + a3*P*v3 + a6*Q*v6 + a8*P*v8 - v6*(a2*v2+a3*v3+a6*v6+a8*v8)/(v2+v3+v6+v8)
a2*Q*v2 + a4*P*v4 + a7*Q*v7 + a8*P*v8 - v7*(a2*v2+a4*v4+a7*v7+a8*v8)/(v2+v4+v7+v8)
a5*Q*v5 + a6*P*v6 + a7*P*v7 + a8*Q*v8 - v8*(a5*v5+a6*v6+a7*v7+a8*v8)/(v5+v6+v7+v8)

где а1, а2, а3, а4, а5, а6, а7, а8, P и Q это коэффициенты.
мой м-файл системы:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
function dx33 = dif33(t,x);
Q=15619375843682503037427100729/15625000000000000000000000000;
P=312406261249325020249757/15625000000000000000000000000;
a1=1;
a2=0.60398;
a3=0.20344;
a4=0.25286;
a5=0.70386;
a6=0.29144;
a7=0.49304;
a8=0.78211;
dx=zeros(8,1);
dx33=[a1*Q*x(1) + a2*P*x(2) + a3*P*x(3) + a4*P*x(4) - x(1)*(a1*x(1)+a2*x(2)+a3*x(3)+a4*x(4))/(x(1)+x(2)+x(3)+x(4));
    a1*P*x(1) + a2*Q*x(2) + a6*P*x(6) + a7*P*x(7) - x(2)*(a1*x(1)+a2*x(2)+a6*x(6)+a7*x(7))/(x(1)+x(2)+x(6)+x(7));
    a1*Q*x(1) + a3*Q*x(3) + a5*P*x(5) + a6*P*x(6) - x(3)*(a1*x(1)+a3*x(3)+a5*x(5)+a6*x(6))/(x(1)+x(3)+x(5)+x(6));
    a1*Q*x(1) + a4*Q*x(4) + a5*P*x(5) + a7*P*x(7) - x(4)*(a1*x(1)+a4*x(4)+a5*x(5)+a7*x(7))/(x(1)+x(4)+x(5)+x(7));
    a3*Q*x(3) + a4*P*x(4) + a5*Q*x(5) + a8*P*x(8) - x(5)*(a3*x(3)+a4*x(4)+a5*x(5)+a8*x(8))/(x(3)+x(4)+x(5)+x(8));
    a2*Q*x(2) + a3*P*x(3) + a6*Q*x(6) + a8*P*x(8) - x(6)*(a2*x(2)+a3*x(3)+a6*x(6)+a8*x(8))/(x(2)+x(3)+x(6)+x(8));
    a2*Q*x(2) + a4*P*x(4) + a7*Q*x(7) + a8*P*x(8) - x(7)*(a2*x(2)+a4*x(4)+a7*x(7)+a8*x(8))/(x(2)+x(4)+x(7)+x(8));
    a5*Q*x(5) + a6*P*x(6) + a7*P*x(7) + a8*Q*x(8) - x(8)*(a5*x(5)+a6*x(6)+a7*x(7)+a8*x(8))/(x(5)+x(6)+x(7)+x(8))];

непосредственно код:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
clear all
clc
 
T=[0 100]; % Интервал интегрирования
        
x0=[10000;0]; % Начальные условия
 
[t,x]=ode15s('dif33',T,x0); %t, x — выходные переменные решателя ode15s
 
plot(t, x(:,1), 'b-', t, x(:,2),'r-',t, x(:,3), 'b-', t, x(:,4),'r-',t, x(:,5), 'b-', t, x(:,6),'r-',t, x(:,7), 'b-', t, x(:,8),'r-')
grid on, 
hold on, 
title('Cистема2'),
legend('V1','V2','V3','V4','V5','V6','V7','V8')
axis([0 1500 -10 15000])

Не могу найти никак ошибки.Если кто-то решал,сталкивался с подобной проблемой,помогите пожалуйста советом!Спасибо за внимание и потраченное время!

Зосима · 21.11.2012, 18:51

Смотри, что этот проказник пишет:
??? Attempted to access x(3); index out of bounds because numel(x)=2.
Это значит, что идет обращение к x(3), в то время как длинна x равна 2!

Действительно, в начальных условиях x0=[10000;0]; - двухэлементный массив, а в функции используется аж до x(8)!
Поэтому могу предложить сделать в качестве начальных условий 8-ми элементный вектор.
Какой именно - сказать не могу, я в дифурах не силен

(сделал x0=[10000;0;0;0;0;0;0;0]; но получилась ерунда и графиков не было)

@Clarisse · 22.11.2012, 09:24 **[ТС]**

Спасибо!!Логичное замечание!!Не досмотрела))Но опять таки да, всё ещё не работает,зато -1 ошибка))

Зосима · 22.11.2012, 09:37

А чего не работает? Что пишет?

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

@Clarisse · 22.11.2012, 09:54 **[ТС]**

А Вы угадали,условия такие и должны были быть(10000,0,0...0).
О чём то предупреждает:
In ode15s at 674
Warning: Matrix is singular, close to singular or badly scaled.
Results may be inaccurate. RCOND = NaN.

Warning: Failure at t=0.000000e+000. Unable to meet integration
tolerances without reducing the step size below the smallest value
allowed (7.905050e-323) at time t.
> In ode15s at 753

и ничего не рисует,может числа маленькие?)

Добавлено через 4 минуты
Вообще ругается на вот эту матрицу,что она плохо масштабируется и т.д. и ругается на шаг:
"Внимание: Если при Т = 0.000000e +000. Не в состоянии удовлетворить интеграции
Допуски без уменьшения размера шага ниже наименьшего значения
разрешено (7.905050e-323) в момент времени Т."
Я как раз не знаю что мне с шагом делать.(

Добавлено через 4 минуты
С рандомной матрицей начальных значений: x0=[1000;1;3;7;1;4;9;1] всё рисуется,но к сожалению мне нужны нули при i=2..8

Зосима · 22.11.2012, 10:17

Если заглянуть во внутрь ф-ции dif33, то можно заметить в последних строках:

Code
1
- x(8)*(a5*x(5)+a6*x(6)+a7*x(7)+a8*x(8))/(x(5)+x(6)+x(7)+x(8))

Очевидно, при x(5)...x(8) = 0 знаменатель равен нулю и получатся то, что получается

Пробовал разные решатели ode45, ode23 и т.д. - тоже самое, что не удивительно.
Видать нужно что-то выдумывать со знаменателем.
*И еще заметил, что в функции dif33 нигде не фигурирует t...

@Clarisse · 22.11.2012, 10:24 **[ТС]**

Дельное замечание!Поговорю с науч. руком.Насчёт t-это как x(t) только там в скобках цифры,чтобы указать x1,x2...x8

Зосима · 22.11.2012, 10:50

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

Насчёт t-это как x(t) только там в скобках цифры,чтобы указать x1,x2...x8

Ясно. Ну как что-то прояснится - говори!

@Clarisse · 22.11.2012, 17:36 **[ТС]**

Заглянула я в документы по курсовой и обнаружила, что в качестве нач. условий можно также использовать вариант:
Vi(0)=N/8,где N= 10000,получается 0 мы больше не имеем и,аллилуйя,оно рисует график)написала руководителю,посмотрю что она скажет)

Зосима · 22.11.2012, 17:40

Молодчика! Поздравляю!

@Clarisse · 22.11.2012, 17:42 **[ТС]**

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

Добавлено через 1 минуту
Это Вам спасибо огромное,ошибку мою заметили,ошибку при вычислении подсказали,без Вас сидела бы голову ломала.Спасибо большое!)

Зосима · 22.11.2012, 18:29

Вот набросал подписки:

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
clear all
clc
 
T = [0 10]; % пределы интегрирования
x0 = 1250 + (1:8)*0; % начальные условия
 
[t,x]=ode15s('dif33',T,x0); %t, x — выходные переменные решателя ode15s
 
% рисуем все
plot(t,x)
 
% подписываем.
% координаты точек подписи каждой кривой.
% подгонялось руками.
k = [10 19 10 11 13 11 12 11];
for i = 1:length(k)
    text(t(k(i)), x(k(i),i), ['\leftarrow V',num2str(i)]);
end
 
hold on
grid on 
title('Cистема2')
legend('V1','V2','V3','V4','V5','V6','V7','V8')
axis([0 10 0 1e4])

Но на мой взгляд, без этих подписей красивее.

Решение системы дифф. уравнений и построение графиков

Зосима · 22.11.2012, 18:31

Кстать, если все-таки решишь не подписывать, а оставить только legend, то вот неплохой вариант отрисовки:

Matlab M
1
2
3
4
5
6
7
8
 plot(t, x(:,1),'b-',...
     t, x(:,2),'k-',...
     t, x(:,3),'g-',...
     t, x(:,4),'r-',...
     t, x(:,5),'c-',...
     t, x(:,6),'m-',...
     t, x(:,7),'y-',...
     t, x(:,8),'m-.')

@Clarisse · 22.11.2012, 18:38 **[ТС]**

Спасибо большое!Первый вариант очень хитрый,не везде понятен лично мне,но спасибо,сохраню,как вариант.
Второй способ сейчас же опробую)

Добавлено через 3 минуты
Наверно без подписей правда аккуратнее.)Зачем эти многоточия после названия цветов?

Зосима · 22.11.2012, 18:39

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

не везде понятен лично мне

Дык спрашивай, еще 20 минут до конца рабочего дня

Добавлено через 24 секунды

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

Зачем эти многоточия после названия цветов?

Перенос на новую строку.

Добавлено через 32 секунды

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

без подписей правда аккуратнее

А без сетки - так хоть на стену вешай!

@Clarisse · 22.11.2012, 18:43 **[ТС]**

отлично)перенос чисто для удобства кода?)
вопросы:
x0 = 1250 + (1:8)*0;сама строка ясна,зачем именно так?
k = [10 19 10 11 13 11 12 11];
for i = 1:length(k)
text(t(k(i)), x(k(i),i), ['\leftarrow V',num2str(i)]); эти строки не совсем поняла,откуда такие числа в массиве и запись
end

axis([0 10 0 1e4]) что за le4 тоже не поняла)

Зосима · 22.11.2012, 18:56

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

перенос чисто для удобства кода?)

Совершенно верно!

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

x0 = 1250 + (1:8)*0;сама строка ясна,зачем именно так?

чтоб не писать вручную восемь раз число 1250)))
Eще можно было: x0 = 1250*ones(8,1);

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

for i = 1:length(k)
text(t(k(i)), x(k(i),i), ['\leftarrow V',num2str(i)]); эти строки не совсем поняла,откуда такие числа в массиве и запись
end

То, что в k - это номера точек времени и функции, как здесь i, только там i было одно, а здесь у нас набор.
Конкретные значения я подбирал руками, чтобы они были на графике и не заслоняли один другой.

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

axis([0 10 0 1e4]) что за le4 тоже не поняла)

1e4=10000

так короче.

@Clarisse · 22.11.2012, 19:01 **[ТС]**

Вроде всё понятно,спасибо за объяснения,просто не привычно)я в матлабе работаю только для курсовой и мало о нём знаю)Сейчас разобралась,спасибо ещё раз огромное за помощь!)

Зосима · 22.11.2012, 19:16

На здоровье!

приятно помочь умной и красивой девушке

Кстать, благодаря Вам я таки немножечко разобрался, как в матлабе решать эти дифуры!

Пасиба!

@Clarisse · 22.11.2012, 20:26 **[ТС]**

Рада что Вы и для себя что-то узнали,помогая мне)))Спасибо что не прошли мимо и потратили на меня своё время!))

Новые блоги и статьи Все статьи Все блоги /
SDL3 для Web (WebAssembly): Идентификация объектов на Box2D v3 - использование userData и событий коллизий 8Observer8 02.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-collision-events-sdl3-c. zip https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11680&d=1772460536 Одним из. . .	Реалии Hrethgir 01.03.2026 Нет, я не закончил до сих пор симулятор. Эта задача сложнее. Не получилось уйти в плавсостав, но оно и к лучшему, возможно. Точнее получалось - но сварщиком в палубную команду, а это значит, в моём. . .	Ритм жизни kumehtar 27.02.2026 Иногда приходится жить в ритме, где дел становится всё больше, а вовлечения в происходящее — всё меньше. Плотный график не даёт вниманию закрепиться ни на одном событии. Утро начинается с быстрых,. . .	SDL3 для Web (WebAssembly): Сборка библиотек: SDL3, Box2D, FreeType, SDL3_ttf, SDL3_mixer и SDL3_image из исходников с помощью CMake и Emscripten 8Observer8 27.02.2026 Недавно вышла версия 3. 4. 2 библиотеки SDL3. На странице официальной релиза доступны исходники, готовые DLL (для x86, x64, arm64), а также библиотеки для разработки под Android, MinGW и Visual Studio. . . .
SDL3 для Web (WebAssembly): Реализация движения на Box2D v3 - трение и коллизии с повёрнутыми стенами 8Observer8 20.02.2026 Содержание блога Box2D позволяет легко создать главного героя, который не проходит сквозь стены и перемещается с заданным трением о препятствия, которые можно располагать под углом, как верхнее. . .	Конвертировать закладки radiotray-ng в m3u-плейлист damix 19.02.2026 Это можно сделать скриптом для PowerShell. Использование . \СonvertRadiotrayToM3U. ps1 <path_to_bookmarks. json> Рядом с файлом bookmarks. json появится файл bookmarks. m3u с результатом. # Check if. . .	Семь CDC на одном интерфейсе: 5 U[S]ARTов, 1 CAN и 1 SSI Eddy_Em 18.02.2026 Постепенно допиливаю свою "многоинтерфейсную плату". Выглядит вот так: https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11617&stc=1&d=1771445347 Основана на STM32F303RBT6. На борту пять. . .	Камера Toupcam IUA500KMA Eddy_Em 12.02.2026 Т. к. у всяких "хикроботов" слишком уж мелкий пиксель, для подсмотра в ESPriF они вообще плохо годятся: уже 14 величину можно рассмотреть еле-еле лишь на экспозициях под 3 секунды (а то и больше),. . .

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	21.11.2012, 18:51
	Смотри, что этот проказник пишет: ??? Attempted to access x(3); index out of bounds because numel(x)=2. Это значит, что идет обращение к x(3), в то время как длинна x равна 2! Действительно, в начальных условиях x0=[10000;0]; - двухэлементный массив, а в функции используется аж до x(8)! Поэтому могу предложить сделать в качестве начальных условий 8-ми элементный вектор. Какой именно - сказать не могу, я в дифурах не силен (сделал x0=[10000;0;0;0;0;0;0;0]; но получилась ерунда и графиков не было) 2

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	22.11.2012, 09:24 [ТС]
	Спасибо!!Логичное замечание!!Не досмотрела))Но опять таки да, всё ещё не работает,зато -1 ошибка)) 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	22.11.2012, 09:37
	А чего не работает? Что пишет? *вводил случайный, без нулей, восьмиэлементный вектор x0 и что-то посчиталось, нарисовалось... 1

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	22.11.2012, 09:54 [ТС]
	А Вы угадали,условия такие и должны были быть(10000,0,0...0). О чём то предупреждает: In ode15s at 674 Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. Warning: Failure at t=0.000000e+000. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t. > In ode15s at 753 и ничего не рисует,может числа маленькие?) Добавлено через 4 минуты Вообще ругается на вот эту матрицу,что она плохо масштабируется и т.д. и ругается на шаг: "Внимание: Если при Т = 0.000000e +000. Не в состоянии удовлетворить интеграции Допуски без уменьшения размера шага ниже наименьшего значения разрешено (7.905050e-323) в момент времени Т." Я как раз не знаю что мне с шагом делать.( Добавлено через 4 минуты С рандомной матрицей начальных значений: x0=[1000;1;3;7;1;4;9;1] всё рисуется,но к сожалению мне нужны нули при i=2..8 0

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	22.11.2012, 10:24 [ТС]
	Дельное замечание!Поговорю с науч. руком.Насчёт t-это как x(t) только там в скобках цифры,чтобы указать x1,x2...x8 0

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	22.11.2012, 17:36 [ТС]
	Заглянула я в документы по курсовой и обнаружила, что в качестве нач. условий можно также использовать вариант: Vi(0)=N/8,где N= 10000,получается 0 мы больше не имеем и,аллилуйя,оно рисует график)написала руководителю,посмотрю что она скажет) Миниатюры 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	22.11.2012, 17:40
	Молодчика! Поздравляю! 1

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	22.11.2012, 17:42 [ТС]
	и вот для этих как раз линий я и хочу подписать их функции,потому что как минимум один цвет совпадает,пунктиры и т.д. не хотелось бы,две линии сливаются и вообще так было бы нагляднее) Добавлено через 1 минуту Это Вам спасибо огромное,ошибку мою заметили,ошибку при вычислении подсказали,без Вас сидела бы голову ломала.Спасибо большое!) 0

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	22.11.2012, 18:38 [ТС]
	Спасибо большое!Первый вариант очень хитрый,не везде понятен лично мне,но спасибо,сохраню,как вариант. Второй способ сейчас же опробую) Добавлено через 3 минуты Наверно без подписей правда аккуратнее.)Зачем эти многоточия после названия цветов? 0

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	22.11.2012, 18:43 [ТС]
	отлично)перенос чисто для удобства кода?) вопросы: x0 = 1250 + (1:8)*0;сама строка ясна,зачем именно так? k = [10 19 10 11 13 11 12 11]; for i = 1:length(k) text(t(k(i)), x(k(i),i), ['\leftarrow V',num2str(i)]); эти строки не совсем поняла,откуда такие числа в массиве и запись end axis([0 10 0 1e4]) что за le4 тоже не поняла) 0

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	22.11.2012, 19:01 [ТС]
	Вроде всё понятно,спасибо за объяснения,просто не привычно)я в матлабе работаю только для курсовой и мало о нём знаю)Сейчас разобралась,спасибо ещё раз огромное за помощь!) 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	22.11.2012, 19:16
	На здоровье! приятно помочь умной и красивой девушке Кстать, благодаря Вам я таки немножечко разобрался, как в матлабе решать эти дифуры! Пасиба! 1

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	22.11.2012, 20:26 [ТС]
	Рада что Вы и для себя что-то узнали,помогая мне)))Спасибо что не прошли мимо и потратили на меня своё время!)) 0