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

勇往直前11

不知道要怎么把ode45求得的数列带入下一个ode45进行计算
: M, \6 v6 t/ j) N2 G% w% z主程序如下：
' Z* ~! g- q8 m%参数初始化
! z# T1 j) [& |8 |clear;
global A B T Q R TH x0;
be=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;
A=[-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];
! h, b  T1 y& w7 M3 UQ=[0 0 0;0 q 0;0 0 0];
) m+ q2 A5 P! |* x$ s, I5 ~R=r;
1 ^+ }9 j3 r6 h0 g; g& b0 h3 hTH=[A -B*B'*inv(R);-Q -A'];
x0=[x10;x10-0;0];
%通过符号运算将相关公式化简方便后面的数值计算
6 H  @. U; y" _, K9 Gsyms s t p11 p12 p13 p22 p23 p33 m1 m2 m3 h1 h2 h3 u1 u2 u3 f1 k1;
/ X0 Y3 p2 R5 I1 _$ Mp=[p11 p12 p13;p12 p22 p23;p13 p23 p33];
  z, E$ K! d0 |/ f2 P( c( y1 B' f4 sp_dot=-A'*p-p*A+p*B*B'*p*inv(R)-Q;
m=[m1;m2;m3];
5 T$ v# f" ~! F7 `m_dot=p*B*B'*m*inv(R)-A'*m;
+ Z1 _+ w! ]# l6 R& @% zh=[h1;h2;h3];
h_dot=p*B*B'*h*inv(R)+p*T-A'*h;
- N* Z& U4 _* r5 f0 iu=[u1 u2 u3];
u_dot=u*B*inv(R)*B'*p-u*A;
f=f1;
f_dot=-u*B*B'*m*inv(R);
5 F. V* S; N  z2 C0 I0 ^+ @k=k1;
: T* }: w7 B; V+ a/ \3 C4 a+ Q3 x  ek_dot=-u*B*B'*h*inv(R)-u*T;
%计算P(t）
tspan=[0,-tf];
0 I3 M% Q$ ^$ q# Z( a: i3 \p0=[0;0;0;0;0;0];
[tp,P]=ode45('DpDt',tspan,p0);
# ]% S9 j. ~6 W0 L/ e" E* Mfigure(1);
. W; n4 Z, [$ d9 |- i/ Bplot(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);
3 H* q- J9 a+ N; r: F3 B%计算M(t）
; i: |! d; O1 \2 a6 Jm0=[0;1;0];
[tm,M]=ode45(@(t,m) DmDt(t,m,tp,P),tspan,m0);
) N% S1 j) ?1 w: Xfigure(2);
plot(tm+tf,M(:,1),'r',tm+tf,M(:,2),'g',tm+tf,M(:,3),'b');
7 n1 {, U6 m7 ~; N* T4 Yylabel('M(t)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
/ O( e' s8 g+ Zset(gca,'XTick',0:0.05:tf,'FontSize',14);
% J7 |9 j1 T% M! \7 H0 w2 s& M$ ]%计算H(t）
h0=[0;0;0];
[th,H]=ode45(@(t,h) DhDt(t,h,tp,P),tspan,h0);
figure(3);
' H: U8 Z# e" K" V3 Bplot(th+tf,H(:,1),'r',th+tf,H(:,2),'g',th+tf,H(:,3),'b');
/ M2 q$ G/ \: R9 ?- a+ P0 f  fylabel('H(t)','FontSize',14);
. k5 ~' h: n, f- J) Axlabel('整个起步过程经历的时间 (s)','FontSize',14);
set(gca,'XTick',0:0.05:tf,'FontSize',14);
%计算U(t）
- T' n( ]  P$ t, g& R6 au0=[0 1 0];
- G- g" P4 N9 @3 a4 ^[tu,U]=ode45(@(t,u) DuDt(t,u,tp,P),tspan,u0);
- ^) i  O8 d8 @figure(4);
% _1 _0 p1 b) Q+ _plot(tu+tf,U(:,1),'r',tu+tf,U(:,2),'g',tu+tf,U(:,3),'b');
ylabel('U(t)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
set(gca,'XTick',0:0.05:tf,'FontSize',14);
%计算F(t）
f0=0;
[tF,F]=ode45(@(t,f) DfDt(t,f,tu,U,tm,M),tspan,f0);
8 F7 X( Z* y- tfigure(5);
plot(tF+tf,F();
ylabel('F(t)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
8 R0 o% W  {  D0 O: Xset(gca,'XTick',0:0.05:tf,'FontSize',14);
# o2 ?/ p  g7 A0 H2 L8 F. c%计算K(t）
2 p+ m9 t# e; Ck0=0;
" f- S. k8 D! t3 a' R7 C4 D[tk,K]=ode45(@(t,k) DkDt(t,f,tu,U,th,H),tspan,k0);
( O: q/ A4 R# Q. |& C: Sfigure(6);
! R% M+ m7 G$ p* M+ Dplot(tk+tf,K();
ylabel('K(t)','FontSize',14);
) E; |4 k; n& Axlabel('整个起步过程经历的时间 (s)','FontSize',14);
$ U# `' q5 u& N5 A6 y" A' \% C9 tset(gca,'XTick',0:0.05:tf,'FontSize',14);
%计算v
$ V% F; M) Q9 m$ ?; AF0=interp1(tF,F(,-tf);
U01=interp1(tu,U(:,1),-tf);
  \. \- R" R$ j& B# ]5 x/ O9 @7 QU02=interp1(tu,U(:,2),-tf);
U03=interp1(tu,U(:,3),-tf);
4 k* j. e0 A0 c9 NK0=interp1(tk,K(:),-tf);
U0=[U01 U02 U03];
2 I. U4 k- y9 _5 v9 ^* ev=inv(F0)*(0-U0*x0-K0);
9 `7 Y0 u% g* j* A  S* c- zfprintf('由公式(26)计算得出的v为：%f\n',v);
- K2 y  b, w9 k) C2 {; e8 O%计算x(t）
  Q( }# N8 M+ Q+ Kx0=[x10;x10-0;0];
tspan=[0,tf];
2 n  F9 S8 F! z# ^& M( i0 i8 n[tx,X]=ode45(@(t,x) DxDt(t,x,tp,P,tm,M,th,H,v,A,B,R,T,tf),tspan,x0);
figure(7);
! S2 {9 O0 m! j# l5 kplot(tx,X(:,2));
ylabel('发动机和离合器从动盘转速差x2(t)  (rad/s)','FontSize',14);
$ y2 U! ?  \# o8 c) Rxlabel('整个起步过程经历的时间 (s)','FontSize',14);
( F! L. Y5 r1 a/ Gset(gca,'XTick',0:0.05:tf,'YTick',0:30:150,'FontSize',14);
figure(8);
+ }1 I" O5 ]' kplot(tx,X(:,3));
4 M$ J( ^; K0 Y( t$ qylabel('离合器正压力x3(t)  (N)','FontSize',14);
xlabel('整个起步过程经历的时间 (s）','FontSize',14);
set(gca,'XTick',0:0.05:tf,'YTick',0:300:1500,'FontSize',14);
4 u2 E. m) u, k0 K7 ?, ?! H! bfigure(9);
# i( a: ]1 y: D: `) K  B* l. }  _plot(tx,X(:,1),'r',tx,X(:,1)-X(:,2),'b');
; \$ I9 F8 Q5 p0 Q- H* Zylabel('发动机和离合器从动盘转速 (rad/s)','FontSize',14);
# U1 i( r* Y! H% p0 exlabel('整个起步过程经历的时间 (s)','FontSize',14);
text(0.5,110,'ωe','FontSize',14);
text(0.5,40,'ωv','FontSize',14);
% x/ L2 O/ _: j5 S0 \6 ?# sset(gca,'XTick',0:0.05:tf,'YTick',0:30:150,'FontSize',14);
& L" E3 Y6 U! K: c3 S% \7 {
, J6 m( r: @5 I
求解P（t）的子程序如下：
function [ p_dot ] =DpDt( t,P )
global A B Q R ;
+ ]/ m: `- H6 o: j2 pp=[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;
p_dot=[pdot(1,1);pdot(1,2);pdot(1,3);pdot(2,2);pdot(2,3);pdot(3,3)];
end

* Q: w: G# h. i" d6 s9 }
求解M（t）的子程序如下：
4 d- S9 e1 C% J0 o; `3 rfunction [ m_dot ] = DmDt( t,m,tp,P )
global A B R;
tp=t;
8 C- {' d# r2 l& ]6 ?m_dot=p*B*inv(R)*B'*m-A'*m;
/ A/ p( ~0 I8 E' L& f* @end
* X! k) y" q( t! r" S$ R

2 M0 m7 D, t* c4 b9 W( O' K结果：
未定义函数或变量 'p'。
出错 DmDt (line 6)
m_dot=p*B*inv(R)*B'*m-A'*m;
5 |7 m2 M# q, Q8 j4 P$ C& o出错 liheqi>@(t,m)DmDt(t,m,tp,P)
出错 odearguments (line 87)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
9 r9 ]+ g* u! d出错 ode45 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn,options, threshold, rtol, normcontrol, normy, hmax, htry,
8 R) x- c9 H, I+ b. Bhtspan, dataType] = odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
; H& k) h& P* H- b# o2 i7 Q  Q6 g2 \出错 liheqi (line 46)
1 X2 `3 _9 J$ G# [+ _9 Y9 _3 X5 c[tm,M]=ode45(@(t,m) DmDt(t,m,tp,P),tspan,m0);
! P- C8 E* ]- |, Q* d# J- F

小小鲁班

帮你顶一下

pTDbn25

楼主，解决了吗？

shuddkk