Info

This question is closed. Reopen it to edit or answer.

Error,Need some help

1 view (last 30 days)
jixie zhuangbei
jixie zhuangbei on 18 Jul 2015
Closed: MATLAB Answer Bot on 20 Aug 2021
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];
%plot(t,x,'r*'); hold on;
k = 3; %多项式的阶次
p=polyfit(t,x,k); %得到拟合的多项式系数
%dt = linspace(0,100,100);
%dx = polyval(p,dt);
%plot(dt,dx);
%xlabel('t');
%ylabel('x');
%title('x-t的拟合曲线');
syms t;%syms t Y;
X = [t^3 t^2 t 1]* p'; %得到X-t的函数
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(t,Y)的符号表达式
Z = matlabFunction(z); %将符号表达式转化为函数
t = linspace(0,180,180);
y(i,6) = Z(t); %Dz-t的方程;dz = Z(3,t);y(i,6) = Z(t);
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,y(:,6)),title('y(6)')
xlabel('t');ylabel('y');
% plot(t,y(:,1),'bo',t,y(:,2),'rx',t,y(:,3),'gv',t,y(:,4),'r-');
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 Comment
Walter Roberson
Walter Roberson on 18 Jul 2015
Change your line
y(i,6) = Z(t);
to
Zt = Z(t);
fprintf('iteration i = %d, number of elements in Zt is %d\n', i, numel(Zt));
y(i,6) = Zt;
This will not prevent the error but it will help you track down the circumstances under which the error occurs. Does it occur on the first iteration, or does it occur after it has been running for some time?

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!