Форум программистов, компьютерный форум, киберфорум
Наши страницы
Matlab
Войти
Регистрация
Восстановить пароль
 
 
Рейтинг 4.63/8: Рейтинг темы: голосов - 8, средняя оценка - 4.63
Gudsaf
104 / 15 / 3
Регистрация: 29.11.2010
Сообщений: 335
1

Один из графиков (ток от времени) не строится корректно! (метод Рунге-Кутты)

19.05.2013, 19:48. Просмотров 1430. Ответов 28
Метки нет (Все метки)

Парни всем привет, не могу понять почему график i(t) всегда идёт по нулям? Т.е. ток всегда нулевой с течением времени. У меня аж 6 вариантов (разные значения данных) и во всех такая проблема. Условие под спойлером.
Кликните здесь для просмотра всего текста

файл-функция системы:
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
function F = dsys(t,I)
% I(1) -> i(t)
% I(2) -> di/dt
% F(1) -> di/dt
% F(2) -> d2i/dt2
 
R0 = 2;
L = 1e-6;
C = 0.001e-6;
f = 1e6;
k = 8e10;
Um = 1;
 
w = 2*pi*f;
R = R0*(1-k*I(1)^2);
 
F = zeros(2,1);
F(1) = I(1);
F(2) = (Um*w*cos(w*t) - 1/C*I(1) - R*I(2) )/L;
программка расчетов:
Matlab M
1
2
3
4
5
6
7
clear, clc
T = [0 5e-6]; 
i0 = [0 0];
[t, I] = ode45('dsys',T,i0);
plot(t,I)
grid on
legend('i(t)','di/dt')
0
Надоела реклама? Зарегистрируйтесь и она исчезнет полностью.
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
19.05.2013, 19:48
Ответы с готовыми решениями:

метод Рунге-Кутты
Методом Рунге-Кутты четвертого порядка на отрезке найти с точностью 10^-5 ...

Метод Рунге-Кутты для решения системы ОДУ
Не могу преобразовать систему в уравнения первого порядка. Помогите пожалуйста...

Метод Рунге-Кутты для решения системы ОДУ
4 дня сижу и пытаюсь понять этот метод(

ОДУ Рунге-Кутты
function =two(f,y0,a,b,eps) % f - имя функции системы % y0 - начальные...

ОДУ методом Рунге-Кутты
y"-4y`+ 5y = 3x y(0)=1.48 y`(0)=3.6 Точн. значение...

28
Зосима
20.05.2013, 11:33
  #2

Не по теме:

*вообще было бы любопытно промоделировать схему в каком-нить MicroCAP :scratch: правда там нелинейное сопротивление...

0
Gudsaf
104 / 15 / 3
Регистрация: 29.11.2010
Сообщений: 335
20.05.2013, 15:50  [ТС] 3
Цитата Сообщение от Зосима Посмотреть сообщение

Не по теме:

*вообще было бы любопытно промоделировать схему в каком-нить MicroCAP :scratch: правда там нелинейное сопротивление...

Не по теме:

я читал в ящике сообщение-оповещение и боялся сейчас увидеть расстрел юзера:)))

0
Зосима
20.05.2013, 16:04
  #4

Не по теме:

Gudsaf, расстрел юзера? та ну, я ж не настолько суровый :-[

Кликните здесь для просмотра всего текста
Один из графиков (ток от времени) не строится корректно! (метод Рунге-Кутты)

0
Gudsaf
104 / 15 / 3
Регистрация: 29.11.2010
Сообщений: 335
20.05.2013, 16:07  [ТС] 5
Цитата Сообщение от Зосима Посмотреть сообщение

Не по теме:

Gudsaf, расстрел юзера? та ну, я ж не настолько суровый :-[

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

Не по теме:

+ 1



А если по теме, то я жду спасителя на белом коне!
0
Зосима
4929 / 3303 / 312
Регистрация: 02.04.2012
Сообщений: 6,207
Записей в блоге: 15
Завершенные тесты: 1
20.05.2013, 17:26 6
Ну, я конечно не спаситель и белого коня у мне неть, но кое-что таки накопал!
Прочитал таки статью, что я тебе кидал по ДУ высших порядков и нашел ошибку:
В системе нужно было прописать F(1) = I(2); ! (это же очевидно по обозначениям di/dt = di/dt)
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
function F = dsys(t,I)
% I(1) -> i(t)
% I(2) -> di/dt
% F(1) -> di/dt
% F(2) -> d2i/dt2
 
R0 = 2;
L = 1e-6;
C = 0.001e-6;
f = 1e6;
k = 0; %8e10;
Um = 1;
 
w = 2*pi*f;
R = R0*(1-k*I(1)^2);
 
F = zeros(2,1);
F(1) = I(2);
F(2) = ( Um*w*cos(w*t) - (1/C)*I(1) - R*I(2) )/L;
Чуток подправил пределы времени:
Matlab M
1
2
3
4
5
6
7
clear, clc
T = [0 6e-6]; 
i0 = [0 0];
[t, I] = ode45('dsys',T,i0);
plot(t,I(:,1))
grid on
legend('i(t)')
И график:
Один из графиков (ток от времени) не строится корректно! (метод Рунге-Кутты)



Время установления режима ~3*10-6c
1
Gudsaf
104 / 15 / 3
Регистрация: 29.11.2010
Сообщений: 335
20.05.2013, 17:34  [ТС] 7
А я пока типарь пишу, спасибо, человеческое такое, знаешь
0
Зосима
4929 / 3303 / 312
Регистрация: 02.04.2012
Сообщений: 6,207
Записей в блоге: 15
Завершенные тесты: 1
20.05.2013, 17:58 8
На здоровье! сам рад, что разобрался.
0
Gudsaf
104 / 15 / 3
Регистрация: 29.11.2010
Сообщений: 335
21.05.2013, 20:12  [ТС] 9
Цитата Сообщение от Зосима Посмотреть сообщение
На здоровье! сам рад, что разобрался.
Зосима, прив, ты тут? Слушай встал я снова с проблемой одной - хочу создать функцию которая принимает множество элементов. На пример так:
Matlab M
1
function F = qweq(t,I,R0,L,C,f,k) %.... далее расчёты
А вызывается она таким образом:
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
        R0 = [2 3 5 3 7.5 1];
        L  = [1e-6 5e-6 10e-6 50e-6 2e-6 10e-6];
        C  = [0.001e-6 0.01e-6 0.1e-6 0.047e-6 0.0047e-6 0.068e-6];
        f  = [1e6 5e6 1e5 2e5 2e6 6e4];
        k  = [8e10 2e14 1e14 5e10 4e15 7e12];
T = [0 6e-6]; 
i0 = [0 0];
i = 1;
while i~=7
    [t, I] = ode45('qweq',T,i0,R0(1,i),L(1,i),C(1,i),f(1,i),k(1,i));
    % ну и там рассчёты ниже
    i = i+1;
end
Видимо я как-то не правильно определяю функцию - он ругает меня.... Я думал там как по подобию Си..

Добавлено через 8 минут
Так ну первую я ошибку нашёл))
Matlab M
1
[t, I] = ode45('qweq',T,i0,R0(1,i),L(1,i),C(1,i),f(1,i),k(1,i));
R0(1,i),L(1,i),C(1,i),f(1,i),k(1,i) - не есть аргументы оде45)

Добавлено через 9 минут
Я надеюсь ты уловил мою идею? я хочу для каждого столбца таблицы (что в задании, оное прикреплено выше) - строить свой график, по сему хочу каждый раз в qweq передавать нужные параметры сопротивления, ёмкости и так далее...

Добавлено через 13 минут
Те вопрос стоит в том, чтобы вызывать функцию qweq каждый раз с разными параметрами... а в оде45, возможно указать в качестве параметров только адрес функции и зарезервированные переменные... Пока кроме как создать 6 функций qweq1 qweq2 qweq3... и так далее со своими параметрами ёмкости, сопротивления и тп никаких идей нет.
0
Зосима
4929 / 3303 / 312
Регистрация: 02.04.2012
Сообщений: 6,207
Записей в блоге: 15
Завершенные тесты: 1
22.05.2013, 10:23 10
Тут можно задействовать глобальные переменные, но мне непонятно, в каком виде ты хочешь получить результат? Т.е. если мы скормим по одному значению данных с каждого массива(R,L,f,U,C,k), то в результате получим два вектора (i(t) и di/dt), да и те, при каждом значении параметров могут иметь разную длинну, поэтому собрать их в одну матрицу будет непросто.
Вот как вариант:
Функция
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function F = dsys(t,I)
global M
% M = [R0, L, C, f, k];
R0 = M(1);
L = M(2);
f = M(4);
Um = 1;
C = M(3);
k = M(5);
 
w = 2*pi*f;
R = R0*(1-k*I(1)^2);
 
F = zeros(2,1);
F(1) = I(2);
F(2) = ( Um*w*cos(w*t) - (1/C)*I(1) - R*I(2) )/L;


Программа расчета
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
global M        
R0 = [2 3 5 3 7.5 1];
L  = [1e-6 5e-6 10e-6 50e-6 2e-6 10e-6];
C  = [0.001e-6 0.01e-6 0.1e-6 0.047e-6 0.0047e-6 0.068e-6];
f  = [1e6 5e6 1e5 2e5 2e6 6e4];
k  = [8e10 2e14 1e14 5e10 4e15 7e12];
 
T = [0 6e-6]; 
i0 = [0 0];
 
for i = 1:length(R)
    M = [R0(i), L(i), C(i), f(i), k(i)];
    [t, I] = ode45('qweq',T,i0);
    % что-то там неведомое
end


*но вот с "k-проблемой" (срыв вычислений после определенного момента t - интеграл не сходится, а шаг уменьшать дальше некуда) пока не выходит справиться получается, что чем больше k, тем меньше время срыва, поэтому синус не успевает раскочегариться...
0
Gudsaf
104 / 15 / 3
Регистрация: 29.11.2010
Сообщений: 335
22.05.2013, 17:17  [ТС] 11
Цитата Сообщение от Зосима Посмотреть сообщение
Тут можно задействовать глобальные переменные, но мне непонятно, в каком виде ты хочешь получить результат? Т.е. если мы скормим по одному значению данных с каждого массива(R,L,f,U,C,k), то в результате получим два вектора (i(t) и di/dt), да и те, при каждом значении параметров могут иметь разную длинну, поэтому собрать их в одну матрицу будет непросто.
Вот как вариант:
Функция
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function F = dsys(t,I)
global M
% M = [R0, L, C, f, k];
R0 = M(1);
L = M(2);
f = M(4);
Um = 1;
C = M(3);
k = M(5);
 
w = 2*pi*f;
R = R0*(1-k*I(1)^2);
 
F = zeros(2,1);
F(1) = I(2);
F(2) = ( Um*w*cos(w*t) - (1/C)*I(1) - R*I(2) )/L;


Программа расчета
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
global M        
R0 = [2 3 5 3 7.5 1];
L  = [1e-6 5e-6 10e-6 50e-6 2e-6 10e-6];
C  = [0.001e-6 0.01e-6 0.1e-6 0.047e-6 0.0047e-6 0.068e-6];
f  = [1e6 5e6 1e5 2e5 2e6 6e4];
k  = [8e10 2e14 1e14 5e10 4e15 7e12];
 
T = [0 6e-6]; 
i0 = [0 0];
 
for i = 1:length(R)
    M = [R0(i), L(i), C(i), f(i), k(i)];
    [t, I] = ode45('qweq',T,i0);
    % что-то там неведомое
end


*но вот с "k-проблемой" (срыв вычислений после определенного момента t - интеграл не сходится, а шаг уменьшать дальше некуда) пока не выходит справиться получается, что чем больше k, тем меньше время срыва, поэтому синус не успевает раскочегариться...
Хорошая идея, то что с глобальными переменными. Да вылет синуса я понял, тоже, пали лс своё. Я пытался использовать функции ode23s ode45s - функции для жёстких функций каких-то... не помогло. В итоге мне просто нужно чтобы рисовало что-то на подобие прикреплённого файла.
Один из графиков (ток от времени) не строится корректно! (метод Рунге-Кутты)
Ведь у нас массив с значениями нормальный, там по возрастанию идут соответствующие значения - причём постепенно возрастают и потом с нуля снова: не пойму почему он не хочет рисовать....

Нет он рисует но он как-то непонятно рисует.... при прорисовке большего интервала по времени получается вот так:
Один из графиков (ток от времени) не строится корректно! (метод Рунге-Кутты)

ТЕ он рисует, как положено (как в моём паинте нарисованном, правда из-за масштаба не видно - там просто полоска) но в конце почему-то херачит вон аж на какой ток - раз в 5000000 больше - потому такая палка страшная.
Надо понять почему он так рассчитывает последнее i
0
Зосима
4929 / 3303 / 312
Регистрация: 02.04.2012
Сообщений: 6,207
Записей в блоге: 15
Завершенные тесты: 1
22.05.2013, 17:43 12
Цитата Сообщение от Gudsaf Посмотреть сообщение
Надо понять почему он так рассчитывает последнее i
скорее всего там в бесконечность улетает (интеграл расходится)
Вопрос в другом - как это преодолеть?
0
Gudsaf
104 / 15 / 3
Регистрация: 29.11.2010
Сообщений: 335
22.05.2013, 18:02  [ТС] 13
Цитата Сообщение от Зосима Посмотреть сообщение
скорее всего там в бесконечность улетает (интеграл расходится)
Вопрос в другом - как это преодолеть?
ты можешь поставить допустим интервал времени от нуля до трёх и будет та же хрень - и выйдет что при трёх он тоже расходится! Но когда мы начертили от 0 до 5 то в трёх он не расходился
0
Зосима
4929 / 3303 / 312
Регистрация: 02.04.2012
Сообщений: 6,207
Записей в блоге: 15
Завершенные тесты: 1
22.05.2013, 21:21 14
Не, если поставить верхний предел поставить, например, 3, то конец графика, этот пик, все-равно будет на ~5.6*10^-6
0
Gudsaf
104 / 15 / 3
Регистрация: 29.11.2010
Сообщений: 335
22.05.2013, 21:25  [ТС] 15
Цитата Сообщение от Зосима Посмотреть сообщение
Не, если поставить верхний предел поставить, например, 3, то конец графика, этот пик, все-равно будет на ~5.6*10^-6
А мне вот сдать надо) и чё делать....вообще почему он туда чертит? там этого значения в таблице нет! он его от балды берёт вообще....

есть какие ещё функции кроме plot?
0
Зосима
4929 / 3303 / 312
Регистрация: 02.04.2012
Сообщений: 6,207
Записей в блоге: 15
Завершенные тесты: 1
23.05.2013, 07:49 16
Здесь дело точно не в plot, а в ode. Эта нелинейность R все ломает, почему-то.
Ты можешь принести два варианта: со скачком, когда k=8e10 и нормальный процесс, когда k <= 100 (тогда время срыва будет больше 6e-6 и мы увидим несколько периодов синусоиды)
Еще я пробовал устанавливать в опциях макс. шаг 1e-20, то тогда программа считает полдня!
Чувствую, что для устранения срыва нужно как-то поиграться с опциями, но я в них мало что понимаю и придерживаюсь принципа: "не трожь, чтоб не сломать".

Кроме того, решение можно попробовать получить в символьном виде, с помошью dsolve.

Не по теме:

Ну и на худой конец, применяя подгониан, в листинге можно написать k = 8e10, а график всунуть при k <= 100, но это для плохих парней, вроде меня :D

0
Gudsaf
104 / 15 / 3
Регистрация: 29.11.2010
Сообщений: 335
23.05.2013, 16:46  [ТС] 17
Цитата Сообщение от Зосима Посмотреть сообщение
Здесь дело точно не в plot, а в ode. Эта нелинейность R все ломает, почему-то.
Ты можешь принести два варианта: со скачком, когда k=8e10 и нормальный процесс, когда k <= 100 (тогда время срыва будет больше 6e-6 и мы увидим несколько периодов синусоиды)
Еще я пробовал устанавливать в опциях макс. шаг 1e-20, то тогда программа считает полдня!
Чувствую, что для устранения срыва нужно как-то поиграться с опциями, но я в них мало что понимаю и придерживаюсь принципа: "не трожь, чтоб не сломать".

Кроме того, решение можно попробовать получить в символьном виде, с помошью dsolve.

Не по теме:

Ну и на худой конец, применяя подгониан, в листинге можно написать k = 8e10, а график всунуть при k <= 100, но это для плохих парней, вроде меня :D

У меня корЯ7 эту хрень при к=90 час считает...
0
Зосима
4929 / 3303 / 312
Регистрация: 02.04.2012
Сообщений: 6,207
Записей в блоге: 15
Завершенные тесты: 1
23.05.2013, 16:48 18
А при k = 10 ?
0
Gudsaf
104 / 15 / 3
Регистрация: 29.11.2010
Сообщений: 335
23.05.2013, 16:51  [ТС] 19
вру, нормально считает

Добавлено через 2 минуты
Цитата Сообщение от Зосима Посмотреть сообщение
А при k = 10 ?
Объясни подробнее что у нас не так, ты сказал что дело всё во времени срыва: если время срыва будет больше 6e-6 и мы увидим несколько периодов синусоиды. А оно зависит от шага, который можно менять, где его менять?
0
Зосима
4929 / 3303 / 312
Регистрация: 02.04.2012
Сообщений: 6,207
Записей в блоге: 15
Завершенные тесты: 1
23.05.2013, 16:55 20

Не по теме:

Цитата Сообщение от Зосима Посмотреть сообщение
с каждого массива(R,L,f,U,C,k)
зря я так расположил переменные :D все дело в этом, инфа 146%



Добавлено через 2 минуты
Цитата Сообщение от Gudsaf Посмотреть сообщение
Объясни подробнее что у нас не так, ты сказал что дело всё во времени срыва: если время срыва будет больше 6e-6 и мы увидим несколько периодов синусоиды. А оно зависит от шага, который можно менять, где его менять?
Оно зависит от k! чем оно больше, тем раньше наступает кирдык!
(если вмазаться, то можно даже получить зависимость tср от k )
0
23.05.2013, 16:55
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
23.05.2013, 16:55

Решение методом Рунге-Кутты
Можно ли без преобразования решить такую систему с помощью MatLab?

Решить методом Рунге-Кутты
j*phi&quot;+r*phi'+mgl(1+pm(y)*sin(w*t+E(y)*sin*phi=-mgl*pm(x)*sin(w*t+E(x)*cos*phi)...

Алгоритм решения ДУ методом Рунге-Кутты
Здравствуйте! Помогите написать алгоритм решения для MatLab методом Рунге-Кутты...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
20
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Рейтинг@Mail.ru