Решение дифференциального уравнения методом Рунге-Кутта

@Altstrelok · Регистрация: 15.12.2015

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

Здравствуйте. Помогите мне, пожалуйста, с такой проблемой. Дано дифференциальное уравнение и необходимо решить методом Рунге-Кутты 4-го порядка с фиксированным шагом, при помощи встроенного оператора и в "ручную". Желательно с подробным описанием, т.к. в этом вопросе еще новенький (пробовал сам, не получается), вторую неделю мучаюсь. Заранее спасибо.
Вот листинг:

Matlab M

clear;
clc;
fur_integr; %здесь задается время, и рассчитываются w1 и Gzol. Pvh и Е - постоянные.
dw1g = @(t,w1g) Gzol.*sqrt(abs(Pvh-(E.*log(w1g./w1))));
w1g=zeros(size(t));
%Начальное значение
w1g(1)=70;
h=0.008; %шаг интегрирования, вообще он в "fur_integr", но чтоб не забыть он здесь.
 
for i=1:length(t)-1
    k1=h*dw1g(t(i), w1g(i));
    k2=h*dw1g(t(i)+h/2, w1g(i)+k1/2);
    k3=h*dw1g(t(i)+h/2, w1g(i)+k2/2);
    k4=h*dw1g(t(i)+h, w1g(i)+k3);
    w1g(i+1)=w1g(i)+(k1+2*k2+2*k3+k4)/6;
        
end

Зосима · 21.12.2015, 21:13

а что не получается?

* попробуй вначале вместо fur_integr; прописать eval fur_integr;
** функция должна быть:

Matlab M

dw1g = @(t,w1g) Gzol.*sqrt(abs(Pvh-(E.*log(w1g./w1)))).*sign( Pvh-(E.*log(w1g./w1)) );

*** решением встроенной функцией:

Matlab M

clear;
clc;
eval fur_integr; %здесь задается время, и рассчитываются w1 и Gzol. Pvh и Е - постоянные.
dw1g = @(t,w1g) Gzol.*sqrt(abs(Pvh-(E.*log(w1g./w1)))).*sign( Pvh-(E.*log(w1g./w1)) );
w1g(1)=70;
h=0.008; %шаг интегрирования, вообще он в "fur_integr", но чтоб не забыть он здесь.
for i=1:length(t)-1
    k1=h*dw1g(t(i), w1g(i));
    k2=h*dw1g(t(i)+h/2, w1g(i)+k1/2);
    k3=h*dw1g(t(i)+h/2, w1g(i)+k2/2);
    k4=h*dw1g(t(i)+h, w1g(i)+k3);
    w1g(i+1)=w1g(i)+(k1+2*k2+2*k3+k4)/6;        
end
[T, W] = ode45(dw1g, [min(t), max(t)], w1g(1) );
plot(T, W, t, w1g)
grid on
legend('ode45','вручную')

@Altstrelok · 22.12.2015, 17:44 **[ТС]**

Вот такую ошибку выдает:
In an assignment A(I) = B, the number of elements in B
and I must be the same.

Error in inet (line 12)
w1g(i+1)=w1g(i)+(k1+2*k2+2*k3+k4)/6;

Добавлено через 11 минут
Error using odearguments (line 91)
@(T,W1G)GZOL.*SQRT(ABS(PVH-(E.*LOG(W1G./W1)))).*SIGN(PVH-(E.*LOG(W1G./W1)))
must return a column vector.

Error in ode45 (line 114)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0,
odeArgs, odeFcn, ...

Зосима · 22.12.2015, 22:27

Altstrelok, без содержимого fur_integr не разберешься

@Altstrelok · 23.12.2015, 17:37 **[ТС]**

Matlab M

clear;
clc;
fn=300;
Nn=80;
tp1=Nn/fn;
tp2=tp1+2*Nn/fn;
t = 0:0.008:0.3;
for i = 1:length(t) 
            N(i)=round(fn.*t(i));
end
 
fis0=0;
Hs=2.5;
fis=N*Hsh-fis0;
 
K0=1/3;
Kmos=14;
Xc=fish*K0/Kmos;
Kc=1.1;
G0=30;
fizol=fis*K0*Kc;
Gzol=G0*abs(fizol);
 
 
wsr=70;
A=20;
L1=5;
L=8.7;
wp=wsr-A*L/2;
w1=wp+A*(L1+Xc);
w2=wp+A*(L-L1-Xc);
 
Kyt=0.5;
E=1000;
Pvh=110;
Pvih=2;

@SSC · 24.12.2015, 07:37

fish - не определен

Массивы
N -> fis -> fizol -> Gzol

Gzol используется при вычислении производной

@Altstrelok · 24.12.2015, 10:23 **[ТС]**

Опечатался малость

Matlab M

Hs=2.5;
fis=N*Hs-fis0;
 
K0=1/3;
Kmos=14;
Xc=fis*K0/Kmos;
Kc=1.1;
G0=30;
fizol=fis*K0*Kc;
Gzol=G0*abs(fizol);
 
 
wsr=70;
A=20;
L1=5;
L=8.7;
wp=wsr-A*L/2;
w1=wp+A*(L1+Xc);
w2=wp+A*(L-L1-Xc);
 
Kyt=0.5;
E=1000;
Pvh=110;
Pvih=2;

@SSC · 24.12.2015, 11:17

А теперь, если этот модуль использовать в модуле из сообщений #1 или #2, то будет неопределен массив t.

@Altstrelok · 24.12.2015, 11:40 **[ТС]**

В смысле не определен? Онже в седьмой строке

@SSC · 24.12.2015, 12:09

Так это в #7 сообщении не весь код а только часть кода из #5 сообщения. Вы просто над нами издеваетесь, мы что тогда проверяем и зачем.
Тогда что Вы сделали для устранения проблемы

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

Массивы
N -> fis -> fizol -> Gzol

Gzol у Вас массив и тогда в этом месте

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

eval fur_integr; %здесь задается время, и рассчитываются w1 и Gzol. Pvh и Е - постоянные.
dw1g = @(t,w1g) Gzol.*sqrt(abs(Pvh-(E.*log(w1g./w1)))).*sign( Pvh-(E.*log(w1g./w1)) );
w1g(1)=70;

что делать с Gzol, какой элемент использовать для вычислений?

@Altstrelok · 24.12.2015, 14:55 **[ТС]**

Поясняю. Сообщение №7 содержит часть исправленного кода в сообщении №5, т.к. при написании сюда, допустил опечатки. Следовательно Gzol вычисляется в нормальном порядке. Теперь вернемся к теме, в "fur_integr" каждый элемент спокойно и правильно вычисляются, а вот далее интегрироваться уравнение не хотит, несмотря отсутствие возражений при написании кода (при написании по правой стороне, при ошибке, замечания подсвечиваются красным или оранжевым цветом). В итоге получается, раз все красиво и гладко почему не интегрируется (высвечивается ошибка указанная выше)?

Добавлено через 6 минут
И если вы думаете, что я ничего не делал или не делаю для решения проблемы, то вы глубоко ошибаетесь! Не в том возрасте чтобы "списывать у кого-то домашку".

@SSC · 24.12.2015, 15:22

Когда Gzol определен в другом программном модуле, какая уж тут подсветка ошибок.
Xc тоже массив, поэтому w1 тоже массив.

Matlab M

clear;
clc;
%eval fur_integr; %здесь задается время, и рассчитываются w1 и Gzol. Pvh и Е - постоянные.
t = 0:0.008:0.3;
Gzol=0.2; % Поставил какую попало
E=1000;
Pvh=110;
wsr=70;
A=20;
L1=5;
L=8.7;
wp=wsr-A*L/2;
Xc=1.2;  % Поставил какую попало
w1=wp+A*(L1+Xc);
 
dw1g = @(t,w1g) Gzol.*sqrt(abs(Pvh-(E.*log(w1g./w1)))).*sign( Pvh-(E.*log(w1g./w1)) );
wg1=zeros(length(t),1);
w1g(1)=70;
h=0.008; %шаг интегрирования, вообще он в "fur_integr", но чтоб не забыть он здесь.
for i=1:length(t)-1
    k1=h*dw1g(t(i), w1g(i));
    k2=h*dw1g(t(i)+h/2, w1g(i)+k1/2);
    k3=h*dw1g(t(i)+h/2, w1g(i)+k2/2);
    k4=h*dw1g(t(i)+h, w1g(i)+k3);
    w1g(i+1)=w1g(i)+(k1+2*k2+2*k3+k4)/6;        
end
[T, W] = ode45(dw1g, [min(t), max(t)], w1g(1) );
plot(T, W, t, w1g)
grid on
legend('ode45','вручную')

@Altstrelok · 24.12.2015, 18:04 **[ТС]**

Хм. Получается он отказывался интегрировать чисто из-за того, что параметры Gzol, Xc, w1 являются массивами? Поскольку эти значения изменяются по времени, они будут образовывать массив, поэтому про такой вариант решения я даже как-то и не подумал.

Хотя, опять же, вопрос остается открытым, пускай и в изменённой форме: "Как тогда проинтегрировать, коль эти значения являются массивами(потому как таких значений может быть не 10 и не 20, а то и от 1000 до миллионов)?". И вот второй вопрос, если не добавлять в код строчку

Matlab M

1	[T, W] = ode45(dw1g, [min(t), max(t)], w1g(1) );

, он конечно решает, но график довольно специфический выводит (горизонтальная прямая из точек), чевой его так "колбасит"?
И можете, на будущее, объяснить мне (как чайнику в этом вопросе), почему он отказывается считать с массивными значениями, ведь массивы-то получаются одинаковыми по размеру?
Извиняюсь за тупые вопросы, но не серчайте пожалуйста, сам постигаю эту науку (нематематик и непрограммист).

@SSC · 25.12.2015, 07:44

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

коль эти значения являются массивами

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

Matlab M

function reshenie
clc
clear all
global R L C U1 U2 W1 W2 Alf1 Alf2
R=20;
L=0.35;
C=500e-6;
U1=550; W1=200*pi; Alf1=-60*pi/180;
U2=100; W2=600*pi; Alf2=0;
% Начальные значения напряжения на С и тока в L
UC_t0=0; 
IL_t0=0;
 
options=odeset('RelTol',1e-7, 'AbsTol',[1e-5 1e-5], 'MaxStep', 0.00001);
[t y]=ode45(@fun_1,[0:0.00001:0.1], [UC_t0 IL_t0], options);
figure(1)
subplot (2,1,1); plot(t,y(:,1));
grid on
xlabel('Время t,c')
ylabel('U_c,В')
subplot (2,1,2); plot(t,y(:,2));
grid on
xlabel('Время t,c')
ylabel('I_L,A')
t(2)
end
 
function dy = fun_1(t,y)
% y(1) -> Uc
% y(2) -> IL
global R L C U1 U2 W1 W2 Alf1 Alf2
dy = zeros(2,1);
E1=U1*sin(W1*t+Alf1)+U2*sin(W2*t+Alf2);
dy(1) =  -y(2)/C;
dy(2) = (E1+y(1)-R*y(2))/L;
end

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

необходимо решить методом Рунге-Кутты 4-го порядка с фиксированным шагом, при помощи встроенного оператора и в "ручную".

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

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

но график довольно специфический выводит

производная не зависит от времени, а только от параметра, да еще под логарифмом, диапазон изменения по времени мал, поэтому и прямая.
Поставь 4 строку
t = 0:0.008:30;
и увидишь другую зависимость.

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

почему он отказывается считать с массивными значениями

Если упрощенно описать задачу которую Вы ставите Matlab-у то выглядит это примерно так:
есть две матрицы 1 строка и n элементов, умнож их нематрично друг на друга и в результате получи массив из 1 элемента.

@Altstrelok · 29.03.2016, 17:59 **[ТС]**

Здравствуйте. Простите за долгий ответ, «немного» не до самообразования было.
В общем, сделал вот так, один раз посчитал, при повторной попытке уже выводит туже ошибку и более не считает. Думаю, посчитало первый раз из-за сбоя программы, хотя вам виднее.

Matlab M

function ishod1
clc
clear all
global fn Nmax Nmin Kn tp1 tp2 fish0 Hsh K0 Kmos Kc G0 wsr A L1 L Kyt E Pvh Pvih  
fn=300; Nmax=130; Nmin=-130; Kn=1; tp1=0.3; tp2=0.9; fish0=0; Hsh=3; K0=1/3;
Kmos=14; Kc=1.1; G0=30; wsr=70; A=20; L1=5; L=8.7; Kyt=0.5; E=1000;
Pvh=110; Pvih=2;
 
w1g_t0=80;
 
options=odeset('RelTol',1e-7, 'AbsTol',[1e-5 1e-5], 'MaxStep', 0.00001);
[t y]=ode45(@fun_1,[0:0.000001:2], [w1g_t0], options);
figure(1)
subplot (1,1,1); plot(t,y(:,1));
grid on
t(2)
end
 
function dy = fun_1(t,y)
% y(1) -> w1g
global fn Nmax Nmin Kn tp1 tp2 fish0 Hsh K0 Kmos Kc G0 wsr A L1 L Kyt E Pvh Pvih w1g
for i = 1:length(t)
    if t(i) <=tp1
        N0=0;
        N(i)=N0+round(Kn*fn.*t(i));
        if N(i)>125
            N(i)=125;
        end
    elseif t(i) <=tp2
        N0=Nmax;
        N(i)=N0-round(Kn*fn.*(t(i)-tp1));
        if N(i)<-125
            N(i)=-125;
        elseif N(i)>125
            N(i)=125;
        end
    else
        N0=Nmin;
        N(i)=N0+round(Kn*fn.*(t(i)-tp2));
        if N(i)<-125
            N(i)=-125;
        elseif N(i)>0
            N(i)=0;
        end
    end
end
fish=N*Hsh-fish0;
Xc=fish*K0/Kmos;
fizol=fish*K0*Kc;
Gzol=G0*abs(fizol);
wp=wsr-A*L/2;
w1=wp+A*(L1+Xc);
w2=wp+A*(L-L1-Xc);
dy = zeros(2,1);
dy(1) =  Gzol.*sqrt(abs(Pvh-(E.*log(w1g./w1)))).*sign( Pvh-(E.*log(w1g./w1)) );
 
end

@SSC · 30.03.2016, 08:18

Я чувствую, бесполезно разбираться с ошибками программы, надо понимать физику того что считается.
Этот вывод я делаю на основании того огорода, который наворочен внутри fun_1.
В первом сообщении есть уравнение (оно конечно выгдядит не как дифференциальное, но допустим что Q это dy/dt)), вопросы по нему:
1. Что такое G, его физический смысл?
2. Что такое w1, его физический смысл?
3. Слева от Q стоит 1 или i - это что такое и важно ли оно для решения задачи?
4. Главное, Вы пытаетесь решать одно диффуравнение, или систему диффуравнений?

Пояснения по солверу ode:

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

function dy = fun_1(t,y)

данная функция предназначена для определения величин производных dy всех параметров решаемой системы в момент времени t при текущих значениях параметров y в данный момент времени.
Внутри этой функции нет ни прошлого, ни будущего, есть только сейчас.
Поэтому

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

for i = 1:length(t)
if t(i) <=tp1

не имеет смысла length(t) всегда равно 1.

Со времен FORTRANa всегда соблюдаю жесткое правило COMMON блоки должны быть одинаковы, иначе почти всегда будешь получать ошибку, причем ее трудно найти. В MATLABе это функция global.
У Вас в основной функции и в fun_1 global отличаются, в результате w1g неопределена, поэтому программа не работает.

Какой физический смысл в округлении (round) внутри fun_1?
Это ведет к очень значительному увеличению времени вычислений.

Введя некоторые (возможно неправомерные поправки) и определив w1g=80;
после продолжительных расчетов получил график

@Altstrelok · 09.04.2016, 08:26 **[ТС]**

Вот такая картина, наиболее полно отражающая физический смысл.
1) G - гидравлическая проводимость дросселирующей щели, которая образуется при перемещении золотника;
2) w1 - геометрический объем полости силового гидроцилиндра, который увеличивается или уменьшается, в зависимости в какой из них втекает жидкость;
3) w1g - гидравлический объем (объем жидкости) в соответствующей полости гидроцилиндра;
Поэтому w1, w1g величины изменяемые во времени.
4) Q - расход жидкости поступающей/истекающей в/из полость(и);
5) В Q1, единица выступает в качестве индекса, показывающая какую полость мы рассматриваем, в данном случае это полость нагнетания.
Да, Q - значение изменяющееся во времени, и своего рода, действительно, является скоростью.
6) Вы правы там несколько дифференциальных уравнений необходимо рассмотреть, но мне достаточно и одного (из личного опыта, остальные диф. уравнения рассматриваются для подтверждения выводов, сделанных из результатов решения этого уравнения).
7) Физический смысл в округлении (round) состоит в том, чтобы упростить/усложнить, соответственно "огрубить"/уточнить, перемещение золотника, от которого зависит значение G. В данном случае, я упростил (грубое приближение) для более простого понимания природы действия.
Если что-то еще непонятно по физике, могу рассказать больше.

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

@SSC · 11.04.2016, 07:58

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

В Q1, единица выступает в качестве индекса

Q1 у Вас вообще нет в уравнениях, а судя по Вашей модели это просто необходимо.

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

Ваша поправка неправомерна.

Судя по вашим описаниям, тогда w1g это y(1).

Судя по модели 1 уравнением Вам не обойтись, как минимум необходимо уравнение по второй полости и уравнения движения массы m, тк она воздействует на поршень препядствуя/способствуя перетокам жидкости.

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

Физический смысл в округлении (round) состоит в том

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

я упростил (грубое приближение) для более простого понимания природы действия.

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

@Altstrelok · 09.05.2016, 20:16 **[ТС]**

Заместо уравнения движения массы можно использовать приведенный момент упругости, обозначен буквой Е (знаю, виноват

), он освобождает от дополнительного/ых уравнения/й. Для второй полости уравнения те же, знак только меняется.
Про

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

вы полностью правы, можно убрать round.
Прошу прощения за долгие ответы. Постараюсь отвечать быстрее.

@Altstrelok 0 / 0 / 0 Регистрация: 15.12.2015 Сообщений: 10
	09.05.2016, 20:16 [ТС]	19
	Заместо уравнения движения массы можно использовать приведенный момент упругости, обозначен буквой Е (знаю, виноват), он освобождает от дополнительного/ых уравнения/й. Для второй полости уравнения те же, знак только меняется. Про система очень жесткая, такие ступеньки вызовут как минимум высокочастотные колебания давления в период прохождения ступеньки вы полностью правы, можно убрать round. Прошу прощения за долгие ответы. Постараюсь отвечать быстрее. 0

@Altstrelok 0 / 0 / 0 Регистрация: 15.12.2015 Сообщений: 10
	22.12.2015, 17:44 [ТС]	3
	Вот такую ошибку выдает: In an assignment A(I) = B, the number of elements in B and I must be the same. Error in inet (line 12) w1g(i+1)=w1g(i)+(k1+2k2+2k3+k4)/6; Добавлено через 11 минут Error using odearguments (line 91) @(T,W1G)GZOL.SQRT(ABS(PVH-(E.LOG(W1G./W1)))).SIGN(PVH-(E.LOG(W1G./W1))) must return a column vector. Error in ode45 (line 114) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ... 0

Зосима 5234 / 3562 / 378 Регистрация: 02.04.2012 Сообщений: 6,470 Записей в блоге: 17
	22.12.2015, 22:27	4
	Altstrelok, без содержимого fur_integr не разберешься 0

@SSC $Эксперт по математике/физике$ 3385 / 1908 / 571 Регистрация: 09.04.2015 Сообщений: 5,350
	24.12.2015, 07:37	6
	fish - не определен Массивы N -> fis -> fizol -> Gzol Gzol используется при вычислении производной 0

@SSC $Эксперт по математике/физике$ 3385 / 1908 / 571 Регистрация: 09.04.2015 Сообщений: 5,350
	24.12.2015, 11:17	8
	А теперь, если этот модуль использовать в модуле из сообщений #1 или #2, то будет неопределен массив t. 0

@Altstrelok 0 / 0 / 0 Регистрация: 15.12.2015 Сообщений: 10
	24.12.2015, 11:40 [ТС]	9
	В смысле не определен? Онже в седьмой строке 0

@SSC $Эксперт по математике/физике$ 3385 / 1908 / 571 Регистрация: 09.04.2015 Сообщений: 5,350
	24.12.2015, 12:09	10
	Так это в #7 сообщении не весь код а только часть кода из #5 сообщения. Вы просто над нами издеваетесь, мы что тогда проверяем и зачем. Тогда что Вы сделали для устранения проблемы Сообщение от SSC Массивы N -> fis -> fizol -> Gzol Gzol у Вас массив и тогда в этом месте Сообщение от Зосима eval fur_integr; %здесь задается время, и рассчитываются w1 и Gzol. Pvh и Е - постоянные. dw1g = @(t,w1g) Gzol.sqrt(abs(Pvh-(E.log(w1g./w1)))).sign( Pvh-(E.log(w1g./w1)) ); w1g(1)=70; что делать с Gzol, какой элемент использовать для вычислений? 0

@Altstrelok 0 / 0 / 0 Регистрация: 15.12.2015 Сообщений: 10
	24.12.2015, 14:55 [ТС]	11
	Поясняю. Сообщение №7 содержит часть исправленного кода в сообщении №5, т.к. при написании сюда, допустил опечатки. Следовательно Gzol вычисляется в нормальном порядке. Теперь вернемся к теме, в "fur_integr" каждый элемент спокойно и правильно вычисляются, а вот далее интегрироваться уравнение не хотит, несмотря отсутствие возражений при написании кода (при написании по правой стороне, при ошибке, замечания подсвечиваются красным или оранжевым цветом). В итоге получается, раз все красиво и гладко почему не интегрируется (высвечивается ошибка указанная выше)? Добавлено через 6 минут И если вы думаете, что я ничего не делал или не делаю для решения проблемы, то вы глубоко ошибаетесь! Не в том возрасте чтобы "списывать у кого-то домашку". 0

@SSC $Эксперт по математике/физике$ 3385 / 1908 / 571 Регистрация: 09.04.2015 Сообщений: 5,350
	30.03.2016, 08:18	16
	Я чувствую, бесполезно разбираться с ошибками программы, надо понимать физику того что считается. Этот вывод я делаю на основании того огорода, который наворочен внутри fun_1. В первом сообщении есть уравнение (оно конечно выгдядит не как дифференциальное, но допустим что Q это dy/dt)), вопросы по нему: 1. Что такое G, его физический смысл? 2. Что такое w1, его физический смысл? 3. Слева от Q стоит 1 или i - это что такое и важно ли оно для решения задачи? 4. Главное, Вы пытаетесь решать одно диффуравнение, или систему диффуравнений? Пояснения по солверу ode: Сообщение от Altstrelok function dy = fun_1(t,y) данная функция предназначена для определения величин производных dy всех параметров решаемой системы в момент времени t при текущих значениях параметров y в данный момент времени. Внутри этой функции нет ни прошлого, ни будущего, есть только сейчас. Поэтому Сообщение от Altstrelok for i = 1:length(t) if t(i) <=tp1 не имеет смысла length(t) всегда равно 1. Со времен FORTRANa всегда соблюдаю жесткое правило COMMON блоки должны быть одинаковы, иначе почти всегда будешь получать ошибку, причем ее трудно найти. В MATLABе это функция global. У Вас в основной функции и в fun_1 global отличаются, в результате w1g неопределена, поэтому программа не работает. Какой физический смысл в округлении (round) внутри fun_1? Это ведет к очень значительному увеличению времени вычислений. Введя некоторые (возможно неправомерные поправки) и определив w1g=80; после продолжительных расчетов получил график Миниатюры 0

@Altstrelok 0 / 0 / 0 Регистрация: 15.12.2015 Сообщений: 10
	09.04.2016, 08:26 [ТС]	17
	Вот такая картина, наиболее полно отражающая физический смысл. 1) G - гидравлическая проводимость дросселирующей щели, которая образуется при перемещении золотника; 2) w1 - геометрический объем полости силового гидроцилиндра, который увеличивается или уменьшается, в зависимости в какой из них втекает жидкость; 3) w1g - гидравлический объем (объем жидкости) в соответствующей полости гидроцилиндра; Поэтому w1, w1g величины изменяемые во времени. 4) Q - расход жидкости поступающей/истекающей в/из полость(и); 5) В Q1, единица выступает в качестве индекса, показывающая какую полость мы рассматриваем, в данном случае это полость нагнетания. Да, Q - значение изменяющееся во времени, и своего рода, действительно, является скоростью. 6) Вы правы там несколько дифференциальных уравнений необходимо рассмотреть, но мне достаточно и одного (из личного опыта, остальные диф. уравнения рассматриваются для подтверждения выводов, сделанных из результатов решения этого уравнения). 7) Физический смысл в округлении (round) состоит в том, чтобы упростить/усложнить, соответственно "огрубить"/уточнить, перемещение золотника, от которого зависит значение G. В данном случае, я упростил (грубое приближение) для более простого понимания природы действия. Если что-то еще непонятно по физике, могу рассказать больше. Ваша поправка неправомерна. Для обывателей (для простоты объяснения или определенных расчетов) жидкость несжимаемая, но если вы помните, при изменении температуры, давления материал изменяет свою форму, и жидкость не исключение. Для простоты, представьте любой газ, который вы сжимаете, к примеру - велосипедный насос, заткнув выходное отверстие вы определенное время спокойно его продолжаете сжимать. Изображения 0

@SSC $Эксперт по математике/физике$ 3385 / 1908 / 571 Регистрация: 09.04.2015 Сообщений: 5,350
	11.04.2016, 07:58	18
	Сообщение от Altstrelok В Q1, единица выступает в качестве индекса Q1 у Вас вообще нет в уравнениях, а судя по Вашей модели это просто необходимо. Сообщение от Altstrelok Ваша поправка неправомерна. Судя по вашим описаниям, тогда w1g это y(1). Судя по модели 1 уравнением Вам не обойтись, как минимум необходимо уравнение по второй полости и уравнения движения массы m, тк она воздействует на поршень препядствуя/способствуя перетокам жидкости. Сообщение от Altstrelok Физический смысл в округлении (round) состоит в том Сообщение от Altstrelok я упростил (грубое приближение) для более простого понимания природы действия. Этого я воообще не понимаю. Вы не на калькуляторе считаете, где набор лишних чисел в основном трата времени практически не увеличивающая точность. В Вашем случае Вы получаете вместо гладкой функции ступенчатую, и это Вы называете упрощением понимания природы действия. А Вы оценили как это скажется на результатах. Ваша система очень жесткая, такие ступеньки вызовут как минимум высокочастотные колебания давления в период прохождения ступеньки. Вам бы периоды открытия/закрытия золотниковой щели пройти, а Вы пораждаете дополнительные возмущения. 0