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

勇往直前11

不知道要怎么把ode45求得的数列带入下一个ode45进行计算
* K5 |6 o% }% j; M主程序如下：
%参数初始化
2 s' X+ z0 t6 u* G6 {; |! }clear;
1 o. ?! _/ [* I: p3 w4 \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];
; ]" G' b! q2 K) ]) \: j& m' y" g% s8 yQ=[0 0 0;0 q 0;0 0 0];
R=r;
# m$ s$ I) _9 H  RTH=[A -B*B'*inv(R);-Q -A'];
x0=[x10;x10-0;0];
+ p+ A8 m+ ?  b  J0 ?2 f# e4 Z%通过符号运算将相关公式化简方便后面的数值计算
- s9 q9 X# C9 \- Isyms s t p11 p12 p13 p22 p23 p33 m1 m2 m3 h1 h2 h3 u1 u2 u3 f1 k1;
5 K" C6 I8 `  d1 f" X) K* Bp=[p11 p12 p13;p12 p22 p23;p13 p23 p33];
p_dot=-A'*p-p*A+p*B*B'*p*inv(R)-Q;
m=[m1;m2;m3];
6 ~8 w3 v/ ^5 i- G& E+ Bm_dot=p*B*B'*m*inv(R)-A'*m;
1 u3 L4 j( \; c6 D2 Zh=[h1;h2;h3];
h_dot=p*B*B'*h*inv(R)+p*T-A'*h;
' N( S: @5 W8 I6 J# Bu=[u1 u2 u3];
u_dot=u*B*inv(R)*B'*p-u*A;
+ r$ y% g/ Z; L" j8 G! v( sf=f1;
* U0 Q" T) U5 W, G0 J3 i+ kf_dot=-u*B*B'*m*inv(R);
7 P, v0 |' ^9 O: O/ g# Y& |$ Qk=k1;
; Y; b8 u" T* l) R2 _k_dot=-u*B*B'*h*inv(R)-u*T;
) f) E/ i3 T; Z8 ?1 b0 L0 Y5 W4 f" n%计算P(t）
( o0 `; |+ b: V; e( d4 Dtspan=[0,-tf];
p0=[0;0;0;0;0;0];
! j6 \) _5 @$ Y6 W! J[tp,P]=ode45('DpDt',tspan,p0);
8 x* [( }: N7 u3 j' {figure(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);
7 P3 h9 g: D9 e3 `* oset(gca,'XTick',0:0.05:tf,'FontSize',14);
%计算M(t）
m0=[0;1;0];
[tm,M]=ode45(@(t,m) DmDt(t,m,tp,P),tspan,m0);
: d' k$ D% F0 `" sfigure(2);
plot(tm+tf,M(:,1),'r',tm+tf,M(:,2),'g',tm+tf,M(:,3),'b');
  G1 ^* c- W4 r- O' zylabel('M(t)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
set(gca,'XTick',0:0.05:tf,'FontSize',14);
%计算H(t）
% Z/ a% W" q/ n: Nh0=[0;0;0];
" l& V% W0 L6 p" s* d; e[th,H]=ode45(@(t,h) DhDt(t,h,tp,P),tspan,h0);
' H, x) z) a! r2 O- f1 ^, tfigure(3);
  Z$ |% H* L6 d1 Z! I0 Dplot(th+tf,H(:,1),'r',th+tf,H(:,2),'g',th+tf,H(:,3),'b');
7 B; j% A% g. b3 a. A" bylabel('H(t)','FontSize',14);
; H) Z: Q) @0 w+ y* Lxlabel('整个起步过程经历的时间 (s)','FontSize',14);
set(gca,'XTick',0:0.05:tf,'FontSize',14);
, ]$ l0 w6 G* M$ X%计算U(t）
$ Z4 L8 ^* C% \' a8 T9 Ku0=[0 1 0];
[tu,U]=ode45(@(t,u) DuDt(t,u,tp,P),tspan,u0);
figure(4);
plot(tu+tf,U(:,1),'r',tu+tf,U(:,2),'g',tu+tf,U(:,3),'b');
ylabel('U(t)','FontSize',14);
3 D- h( c8 @; E# s# ^$ F1 W6 Axlabel('整个起步过程经历的时间 (s)','FontSize',14);
set(gca,'XTick',0:0.05:tf,'FontSize',14);
# G, u  {( A" e5 k! e) w8 D%计算F(t）
f0=0;
[tF,F]=ode45(@(t,f) DfDt(t,f,tu,U,tm,M),tspan,f0);
figure(5);
: [% s, R- ]6 B1 Tplot(tF+tf,F();
3 U$ d+ N5 o" O  m: E# t! f" _ylabel('F(t)','FontSize',14);
! ?+ B5 E% X4 X3 ^& b1 kxlabel('整个起步过程经历的时间 (s)','FontSize',14);
0 D5 M. [, ]; z0 Y  iset(gca,'XTick',0:0.05:tf,'FontSize',14);
5 n$ _3 p( ?% ?%计算K(t）
k0=0;
$ b: `8 o9 {8 o& H[tk,K]=ode45(@(t,k) DkDt(t,f,tu,U,th,H),tspan,k0);
figure(6);
plot(tk+tf,K();
" c# G0 o" h+ b1 D8 t! _" x! E4 qylabel('K(t)','FontSize',14);
xlabel('整个起步过程经历的时间 (s)','FontSize',14);
! i+ Q$ s* Q/ d9 x8 Eset(gca,'XTick',0:0.05:tf,'FontSize',14);
%计算v
F0=interp1(tF,F(,-tf);
1 e  o9 C6 c& U3 a4 ^  \/ w+ jU01=interp1(tu,U(:,1),-tf);
8 d7 q- x7 r  T6 @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);
fprintf('由公式(26)计算得出的v为：%f\n',v);
%计算x(t）
3 i  E1 p1 q3 {1 r1 K' o7 m0 ^x0=[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);
plot(tx,X(:,2));
1 L" s6 [' k5 s; W  @) L# sylabel('发动机和离合器从动盘转速差x2(t)  (rad/s)','FontSize',14);
. _; @  }6 e$ a2 [4 R: s' yxlabel('整个起步过程经历的时间 (s)','FontSize',14);
set(gca,'XTick',0:0.05:tf,'YTick',0:30:150,'FontSize',14);
figure(8);
plot(tx,X(:,3));
$ ~. X8 q5 T1 M2 {. T& F6 Z! s1 P4 lylabel('离合器正压力x3(t)  (N)','FontSize',14);
/ H8 [; X3 c: Ixlabel('整个起步过程经历的时间 (s）','FontSize',14);
set(gca,'XTick',0:0.05:tf,'YTick',0:300:1500,'FontSize',14);
figure(9);
plot(tx,X(:,1),'r',tx,X(:,1)-X(:,2),'b');
ylabel('发动机和离合器从动盘转速 (rad/s)','FontSize',14);
0 F3 o1 j( \$ k! p8 e8 }4 wxlabel('整个起步过程经历的时间 (s)','FontSize',14);
text(0.5,110,'ωe','FontSize',14);
text(0.5,40,'ωv','FontSize',14);
- c$ Y) a% ~7 r% K, vset(gca,'XTick',0:0.05:tf,'YTick',0:30:150,'FontSize',14);

* L3 `. f3 g; v9 D4 ^$ D7 l( ?$ x
求解P（t）的子程序如下：
function [ p_dot ] =DpDt( t,P )
global A B Q R ;
p=[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)];
3 s$ T* L$ }2 ^6 opdot=-A'*p-p*A+p*B*B'*p*inv(R)-Q;
1 y) V2 g# s: a! T& O( Jp_dot=[pdot(1,1);pdot(1,2);pdot(1,3);pdot(2,2);pdot(2,3);pdot(3,3)];
end


! `, r& J6 y" _+ Q求解M（t）的子程序如下：
! K7 s1 j6 x7 @function [ m_dot ] = DmDt( t,m,tp,P )
global A B R;
2 d  e# F8 X; _3 o. A7 btp=t;
1 Y" T, |7 f! y, p) Pm_dot=p*B*inv(R)*B'*m-A'*m;
end
8 x7 w& J: b, z" ~- O
. @7 r$ F! E8 M
结果：
未定义函数或变量 'p'。
, i  n* J4 ]& i5 v4 q% P出错 DmDt (line 6)
1 x' V! g6 f1 o- Y# i3 um_dot=p*B*inv(R)*B'*m-A'*m;
$ W) n) g! j5 ^. r. A7 w, x/ e% {: {0 \% Q出错 liheqi>@(t,m)DmDt(t,m,tp,P)
3 s: U* }) q0 {! r: K+ t. G出错 odearguments (line 87)
; i. C+ |, u8 }+ z* P: ^f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
出错 ode45 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn,options, threshold, rtol, normcontrol, normy, hmax, htry,
htspan, dataType] = odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
出错 liheqi (line 46)
! Y$ x2 x" U( g# n$ e[tm,M]=ode45(@(t,m) DmDt(t,m,tp,P),tspan,m0);

小小鲁班

帮你顶一下

pTDbn25

楼主，解决了吗？

shuddkk