Wrong，all the time（Equation（6） is y(i,6)， not dy(6)）

Question

jixie zhuangbei 2015-7-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/230728-wrong-all-the-time-equation-6-is-y-i-6-not-dy-6

评论： jixie zhuangbei 2015-7-23

Equation：

dy(1)=-1.918298553*y(3)*y(4)-121.6697369*y(5)*y(2)+0.006472085*y(2)*y(2)+15.25250926*y(5)*y(5)-0.518363603*sin(2.047*t)+0.001124759;
dy(2)=0.007229182*y(5)*y(1)-0.013867729*y(3)-0.005151943*sin(2.047*t)+33.43424564*y(4)*y(5)-0.092169794*y(3)*y(2)-0.698266828;
dy(3)=72.10986245*y(2)+0.52129529*y(4)*y(1)+0.025471074*y(3)+921.886526*y(4)/y(1)-0.47870471*y(4)+0.025471074*y(1)*t-0.38220722*cos(2.047*t)-4.62279911;
dy(4)=-57.37263009*y(5)+0.001053501*y(3)*y(1)+20.28825016/y(1)-0.001053501*y(3)+0.064741605*y(4)+20.28825016*t-0.006201915*sin(2.047*t)-3651.885374303154;
dy(5)=0.017284293*y(4)+0.00278644*y(1)*y(2)-0.551218454*y(2)*y(5)*y(5)+0.010839281*y(2)*y(2)*y(5)+0.020353114*cos(2.047*t)+0.110594984;
y(i,6) = -178075801.6*y(i,4)*y(i,5)+X*y(i,3)*y(i,2)+1119396.815471592+7481.104004*sin(2.047*t(i));

In order to solve the equations， I edited 2 M files. In equation (6), X is coefficient which varies with t(t is time)

【ODE451_main】

clc;clear
tspan=[0,180];
y0=[1e-5;0;0;0;0];
[t,y]=ode23('ODE45_fun',tspan,y0);
[m,n]=size(y);
for i=1:m
    t=[0 25 50 70 80 100 180];
    x=[0 5  7  9  10 11  14];
k = 3;   
p=polyfit(t,x,k); 
syms t;%syms t Y;
X = [t^3 t^2 t 1]* p'; 
z =-178075801.6*y(i,4)*y(i,5)+X*y(i,3)*y(i,2)+1119396.815471592+7481.104004*sin(2.047*t(i));
Z = matlabFunction(z);  
end
data=[t,y]; 
save ODE45_data.txt data -ascii
subplot(2,3,1),plot(t,y(:,1)),title('y(1)')
xlabel('t');ylabel('y');
subplot(2,3,2),plot(t,y(:,2)),title('y(2)')
xlabel('t');ylabel('y');
subplot(2,3,3),plot(t,y(:,3)),title('y(3)')
xlabel('t');ylabel('y');
subplot(2,3,4),plot(t,y(:,4)),title('y(4)')
xlabel('t');ylabel('y');
subplot(2,3,5),plot(t,y(:,5)),title('y(5)')
xlabel('t');ylabel('y');
subplot(2,3,6),plot(t,Z),title('Z')  
xlabel('t');ylabel('Z');
grid on

【ODE45_fun】

function  dy=ODE45_fun(t,y)
dy(1)=-1.918298553*y(3)*y(4)-121.6697369*y(5)*y(2)+0.006472085*y(2)*y(2)+15.25250926*y(5)*y(5)-0.518363603*sin(2.047*t)+0.001124759;
dy(2)=0.007229182*y(5)*y(1)-0.013867729*y(3)-0.005151943*sin(2.047*t)+33.43424564*y(4)*y(5)-0.092169794*y(3)*y(2)-0.698266828;
dy(3)=72.10986245*y(2)+0.52129529*y(4)*y(1)+0.025471074*y(3)+921.886526*y(4)/y(1)-0.47870471*y(4)+0.025471074*y(1)*t-0.38220722*cos(2.047*t)-4.62279911;
dy(4)=-57.37263009*y(5)+0.001053501*y(3)*y(1)+20.28825016/y(1)-0.001053501*y(3)+0.064741605*y(4)+20.28825016*t-0.006201915*sin(2.047*t)-3651.885374303154;
dy(5)=0.017284293*y(4)+0.00278644*y(1)*y(2)-0.551218454*y(2)*y(5)*y(5)+0.010839281*y(2)*y(2)*y(5)+0.020353114*cos(2.047*t)+0.110594984;
dy=[dy(1);dy(2);dy(3);dy(4);dy(5)];

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Stephen23 2015-7-22

编辑：Stephen23 2015-7-22

Please format your code properly to make it readable. This means doing all of these steps:

edit your question.
remove all empty lines form the code
select the code
click the {} Code button that you will find above the textbox.

Currently your code is unreadable. Also you might like to actually ask a question.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2015-7-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/230728-wrong-all-the-time-equation-6-is-y-i-6-not-dy-6#answer_186831

在 MATLAB Online 中打开

Instead of your for-loop, use

t0=[0 25 50 70 80 100 180];
x=[0 5  7  9  10 11  14];
k = 3;
p=polyfit(t0,x,k);
X=polyval(p,t);
y(:,6)=-178075801.6*y(:,4).*y(:,5)+X.*y(:,3).*y(:,2)+1119396.815471592+7481.104004*sin(2.047*t(:));

Best wishes

Torsten.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Torsten 2015-7-23

Technically, your code works, but ODE23 is not able to integrate your equations properly (only up to t=1e-23 s). You should check them again for errors.

I guess that the division by y(1) in equations 3 and 4 causes the trouble.

Best wishes

Torsten.

jixie zhuangbei 2015-7-23