Решение дифференциальных уравнений

function dy=odefun10(t,y);
N=3;
J=12800;
q=350;
S=38;
l=3.7;
ks=-0.005;
ko=-0.01;
ka=-0.01;
a=3;
kz=-0.08671;
 
%Расчеты
m1=y(1);
omega=y(2);
tetta=y(3);
m=2*(N-tetta);
sigma=(J/(q*S*l*ks))*(3*(m-omega)+m1)-(ko+ka*a+kz*omega)\ks;
 
dy=zeros(3,1);
dy(1)= m;
dy(2)=((q*S*l)\J)*(ko+ka*a+ks*sigma+kz*omega);
dy(3)=omega;
end

[Y,T]=ode45('odefun10', [0 10], [0 0 0]);
plot(Y,T)

function OdeCaller2
N=3;
J=12800;
q=350;
S=38;
l=3.7;
ks=-0.005;
ko=-0.01;
ka=-0.01;
a=3;
kz=-0.08671;
 
[T,Y]=ode45(@odefun10, [0 4], [0 0]);
 
m1=-2*Y(:,1);
plot(T,m1,'r-')
hold on
plot(T,Y)
grid on
legend('m1','\omega','\upsilon')
function dy=odefun10(t,y)
%Расчеты
omega=y(1);
tetta=y(2);
m=2*(N-tetta);
m1=-2*omega;
% Что бы понять откуда я взял такую формулу введите следующие 2 строки в коммандном
% окне:
% syms N tetta(t)
% m1=diff(2*(N-tetta),t)
% а затем взгляните на последнее дифференциальное уравнение на вашей
% картинке
dy=zeros(2,1);
 
sigma=(J/(q*S*l*ks))*(3*(m-omega)+m1)-(ko+ka*a+kz*omega)/ks;
 
dy(1)=((q*S*l)/J)*(ko+ko*a+ks*sigma+kz*omega);
dy(2)=omega;
end
end

for i=1:6;
y(i)=rand;
S=S+y(i);
end
E=(S+3)\6;

function OdeCaller5
N=0;
S=0
 
[T,Y]=ode45(@odefun10, [0 4], [0 0]);
Uk=N;
 
dy(1)=L2;
dy(2)=sigma1;
 
for i=1:6;
Q(i)=2*rand-1;
S=S+Q(i);
end
E=(S+3)\6;
 
%b=exp(-0.1*дельта t); где дельта t-значение шага интегрирования
Vy=b*Uk + 2*sqrt(1-b^2)*E;
 
L1=0.06*Vy;
 
dy=zeros(2,1);
dy(1)=0.01*Vy-5*L2
 
L=L1+L2;
 
dy(2)=L;
 
N=Vy
 
end
end

function dy=odefun10(t,y)

function OdeCaller5
N=0;
S=0;
 
[T,Y]=ode45(@odefun11, [0 10], [0 0]);
 
plot(T,Y(:,2))
 
 
function dy=odefun11(t,y)
Uk=N;
L2=y(1);
 
for i=1:6;
Q(i)=2*rand-1;
S=S+Q(i);
end
E=(S+3)/6;
 
%b=exp(-0.1*дельта t); где дельта t-значение шага интегрирования
 
Vy=b*Uk + 2*sqrt(1-b^2)*E;
 
L1=0.006*Vy;
 
dy=zeros(2,1);
dy(1)=0.01*Vy-5*L2;
 
L=L1+L2;
 
dy(2)=L;
 
N=Vy;
 
end
end

[T,Y]=ode45(@odefun10, 0:Tstep:4, [0 0]);
% где Tstep это шаг интегрирования

function OdeCaller
ik=0; % Инициализация счётчика
...
Tspan=0:Tstep:4;
VyMem=zeros(1,numel(Tspan));% Заготавливаем память под VyMem
[T,Y]=ode45(@odefun10,Tspan , [0 0]);
  function dy=odefun11(t,y)
  ik=ik+1; %Считаем сколько раз вызвана функция odefun11
  ...
  VyMem(ik)=Vy;% Собираем данные в массив
  end
end

function OdeCaller5
ik=0;% Инициализация счётчика
tstep=0.001;
Tspan=0:tstep:40;
N=0;
VyMem=zeros(1,numel(Tspan));% Заготавливаем память под VyMem
[T,Y]=ode45(@odefun11,Tspan, [0 0]);
plot(T,Y)
 
 
function dy=odefun11(t,y)
ik=ik+1; %Считаем сколько раз вызвана функция odefun11
Uk=N;
L2=y(1);
E=randn;
 
b=exp(-0.1*tstep);
Vy=b*Uk + 2*sqrt(1-b^2)*E;
 
L1=0.06*Vy;
 
dy=zeros(2,1);
dy(1)=0.001*Vy-5*L2;
 
L=L1+L2;
 
dy(2)=L;
N=Vy;
VyMem(ik)=Vy;
end
end