ode45输出数列带入下一个ode45进行计算

勇往直前11

不知道要怎么把ode45求得的数列带入下一个ode45进行计算
主程序如下：
9 v7 N: J& [% H. J, D%参数初始化
clear;
9 [9 Y8 F/ G- lglobal A B T Q R TH x0;
  j% V, u: R& E% b; ]4 c$ E8 ^( ube=0.03; bv= 0.03; Ie=0.2; Iv=0.7753; k=0.098; ud=0.4; Fnmax=5000;
Tin=100; Tl=10; tf=0.65; q=1000; r=1; umax=8000; x10=95; x1min=70;
+ G& D6 S2 t5 P' f% zA=[-be/Ie 0 -k/Ie;(-be/Ie+bv/Iv) -bv/Iv -(k/Ie+k/Iv);0 0 0];
B=[0;0;1];
T=[Tin/Ie;(Tin/Ie+Tl/Iv);0];
Q=[0 0 0;0 q 0;0 0 0];
$ g) M" }8 e3 X. H  r) j+ uR=r;
4 k2 ?5 u, |9 M) j* R( ]& ?' bTH=[A -B*B'*inv(R);-Q -A'];
x0=[x10;x10-0;0];
/ g- a' T. @  q! N- f$ V%通过符号运算将相关公式化简方便后面的数值计算
syms s t p11 p12 p13 p22 p23 p33 m1 m2 m3 h1 h2 h3 u1 u2 u3 f1 k1;
' K2 Y: y/ D+ N8 k' }# }$ @p=[p11 p12 p13;p12 p22 p23;p13 p23 p33];
$ U/ u& F+ ~$ \; v* u0 A7 |p_dot=-A'*p-p*A+p*B*B'*p*inv(R)-Q;
3 R4 M. Y* |! T3 h* q$ T8 dm=[m1;m2;m3];
/ {0 D% W0 r( F; v: gm_dot=p*B*B'*m*inv(R)-A'*m;
h=[h1;h2;h3];
. }: J7 ]9 C6 H3 `' [5 wh_dot=p*B*B'*h*inv(R)+p*T-A'*h;
u=[u1 u2 u3];
u_dot=u*B*inv(R)*B'*p-u*A;
f=f1;
f_dot=-u*B*B'*m*inv(R);
9 \: t! |9 G6 i+ ck=k1;
k_dot=-u*B*B'*h*inv(R)-u*T;
%计算P(t）
tspan=[0,-tf];
2 I# \) B1 @" c# Pp0=[0;0;0;0;0;0];
4 a0 [5 ~" E: [4 l: X+ _[tp,P]=ode45('DpDt',tspan,p0);
# c6 c% M% y7 S* hfigure(1);
plot(tp+tf,P(:,1),'r',tp+tf,P(:,2),'g',tp+tf,P(:,3),'b',tp+tf,P(:,4),'c',tp+tf,P(:,5),'m',tp+tf,P(:,6),'k');
ylabel('黎卡提微分方程的解 P(t)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
set(gca,'XTick',0:0.05:tf,'FontSize',14);
%计算M(t）
" k& V7 H1 A( c- [* fm0=[0;1;0];
[tm,M]=ode45(@(t,m) DmDt(t,m,tp,P),tspan,m0);
1 H* m4 |8 ~/ Y! ?figure(2);
! v& `! ?* K# z' F9 wplot(tm+tf,M(:,1),'r',tm+tf,M(:,2),'g',tm+tf,M(:,3),'b');
4 W1 _" Y+ x6 P; U$ }* [ylabel('M(t)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
4 ~; s5 a; H3 c  nset(gca,'XTick',0:0.05:tf,'FontSize',14);
& Q$ @* F! o  f/ O%计算H(t）
h0=[0;0;0];
9 F& B' n# i0 _1 j) a& V, S1 I2 U7 p[th,H]=ode45(@(t,h) DhDt(t,h,tp,P),tspan,h0);
figure(3);
% i, b( U# p4 x' oplot(th+tf,H(:,1),'r',th+tf,H(:,2),'g',th+tf,H(:,3),'b');
* o1 I' y+ G4 p- G9 ^ylabel('H(t)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
set(gca,'XTick',0:0.05:tf,'FontSize',14);
3 [" J' H9 u1 x/ y! {( L%计算U(t）
u0=[0 1 0];
[tu,U]=ode45(@(t,u) DuDt(t,u,tp,P),tspan,u0);
figure(4);
: _. c( R$ B& v" Splot(tu+tf,U(:,1),'r',tu+tf,U(:,2),'g',tu+tf,U(:,3),'b');
- h" S$ g1 q* \" Q) d5 \ylabel('U(t)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
set(gca,'XTick',0:0.05:tf,'FontSize',14);
%计算F(t）
" P9 h0 K! s: ^, f, F, G) w; W5 rf0=0;
[tF,F]=ode45(@(t,f) DfDt(t,f,tu,U,tm,M),tspan,f0);
figure(5);
& J2 T' a) Q1 ], }6 q1 l" S5 Nplot(tF+tf,F();
ylabel('F(t)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
set(gca,'XTick',0:0.05:tf,'FontSize',14);
%计算K(t）
! h8 R% O" P9 ]3 Hk0=0;
[tk,K]=ode45(@(t,k) DkDt(t,f,tu,U,th,H),tspan,k0);
8 Q9 O: `; Q! f- J( n2 Cfigure(6);
! ?5 m. A; L- ~3 c! eplot(tk+tf,K();
ylabel('K(t)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
8 x; }. j' A, C# q6 bset(gca,'XTick',0:0.05:tf,'FontSize',14);
%计算v
F0=interp1(tF,F(,-tf);
U01=interp1(tu,U(:,1),-tf);
U02=interp1(tu,U(:,2),-tf);
U03=interp1(tu,U(:,3),-tf);
K0=interp1(tk,K(:),-tf);
U0=[U01 U02 U03];
v=inv(F0)*(0-U0*x0-K0);
' j7 O3 @  C4 ^+ e# B0 Pfprintf('由公式(26)计算得出的v为：%f\n',v);
%计算x(t）
8 L9 H! l) i6 j% N6 nx0=[x10;x10-0;0];
tspan=[0,tf];
[tx,X]=ode45(@(t,x) DxDt(t,x,tp,P,tm,M,th,H,v,A,B,R,T,tf),tspan,x0);
figure(7);
6 C3 t; d* \8 Dplot(tx,X(:,2));
ylabel('发动机和离合器从动盘转速差x2(t)  (rad/s)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
set(gca,'XTick',0:0.05:tf,'YTick',0:30:150,'FontSize',14);
5 Z& B1 ~+ j: tfigure(8);
, w3 q( B7 w  u8 bplot(tx,X(:,3));
ylabel('离合器正压力x3(t)  (N)','FontSize',14);
3 o8 J& F( w* H) \4 ?9 @xlabel('整个起步过程经历的时间 (s）','FontSize',14);
set(gca,'XTick',0:0.05:tf,'YTick',0:300:1500,'FontSize',14);
figure(9);
/ h7 i3 t0 |% `# R( Aplot(tx,X(:,1),'r',tx,X(:,1)-X(:,2),'b');
9 ^. u! J( H; \6 w- O& B/ wylabel('发动机和离合器从动盘转速 (rad/s)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
4 m' ~* Z" F+ C2 G: Q3 s3 W) e6 Mtext(0.5,110,'ωe','FontSize',14);
text(0.5,40,'ωv','FontSize',14);
& `  Z) d" A; ~" H% x4 Sset(gca,'XTick',0:0.05:tf,'YTick',0:30:150,'FontSize',14);
  y6 {- c  @2 o7 m& V
3 r6 k) [! K  p/ a8 ~' M; x; \; |
求解P（t）的子程序如下：
$ C+ {: T- h& [! @function [ p_dot ] =DpDt( t,P )
, h; ?: g* x, v4 `- s! R8 sglobal A B Q R ;
& k9 a% L, D) s( l. q( |, g7 op=[P(1,1) P(2,1) P(3,1)(2,1) P(4,1) P(5,1)(3,1) P(5,1) P(6,1)];
pdot=-A'*p-p*A+p*B*B'*p*inv(R)-Q;
) n# H. ^- M( T' ?p_dot=[pdot(1,1);pdot(1,2);pdot(1,3);pdot(2,2);pdot(2,3);pdot(3,3)];
& A% X( z- {. X) Q( Lend


8 H0 U. y) N6 Q% Z求解M（t）的子程序如下：
function [ m_dot ] = DmDt( t,m,tp,P )
global A B R;
tp=t;
# [. p! n) X! @: q: Km_dot=p*B*inv(R)*B'*m-A'*m;
2 `, q$ `# X, g5 O: I* u, N0 @# Kend


结果：
; Q' T, W3 E* X& |; n! N未定义函数或变量 'p'。
出错 DmDt (line 6)
! m( u  {! m6 V% b0 W( Tm_dot=p*B*inv(R)*B'*m-A'*m;
出错 liheqi>@(t,m)DmDt(t,m,tp,P)
% p' [5 o/ J# L0 U2 y9 m# _出错 odearguments (line 87)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
1 v3 k$ d' N; w& M6 H: D出错 ode45 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn,options, threshold, rtol, normcontrol, normy, hmax, htry,
+ x$ E4 Z. ?# I$ ]& s) E% ohtspan, dataType] = odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
* |8 j# Q3 D  i6 n( v出错 liheqi (line 46)
[tm,M]=ode45(@(t,m) DmDt(t,m,tp,P),tspan,m0);

小小鲁班

帮你顶一下

pTDbn25

楼主，解决了吗？

shuddkk