我用ode45求解了一个二阶线性非齐次方程,得到x1和x2=dx1/dt与时间t 的关系图,程序如下:
function linearResponse
m0=0.46;r=13.5; w=200*2*3.14;m=3.76+0.46; f=578;k=10;
tspan=linspace(0,8,100); sampint=tspan(2);
options=odeset('RelTol',1e-8,'AbsTol',[1e-8,1e-8]);
[t,x]=ode45(@ForcedOscillatorl,tspan,[0 0]',options,m0,r,w,m,f,k);
figure(1);
subplot(2,1,1); plot(t,x(:,1)); axis([0 8,-2,2]); xlabel('\tau'); ylabel('x(\tau)');
subplot(2,1,2); plot(t,x(:,2)); axis([0 8,-2000,2000]); xlabel('\tau'); ylabel('v(\tau)');
function xdot=ForcedOscillatorl(t,x,m0,r,w,m,f,k)
xdot=[x(2);-f*x(2)/(m0+m)-k*x(1)/(m0+m)+m0*r*cos(w*t)*w^2/(m0+m)]
可是我想还得到dx2/dt与时间t 的关系图该怎么办呢?
用diff(x(,:2))好像是不行的,程序:
function linearResponse
m0=0.46;r=13.5; w=200*2*3.14;m=3.76+0.46; f=578;k=10;
tspan=linspace(0,8,100); sampint=tspan(2);
options=odeset('RelTol',1e-8,'AbsTol',[1e-8,1e-8]);
[t,x]=ode45(@ForcedOscillatorl,tspan,[0 0]',options,m0,r,w,m,f,k);
figure(1);
subplot(3,1,1); plot(t,x(:,1)); axis([0 8,-2,2]); xlabel('\tau'); ylabel('x(\tau)');
subplot(3,1,2); plot(t,x(:,2)); axis([0 8,-2000,2000]); xlabel('\tau'); ylabel('v(\tau)');
x3=diff(x(:,2));
subplot(3,1,3); plot(t,x3); axis([0 8,-2000,2000]); xlabel('\tau'); ylabel('a(\tau)');
function xdot=ForcedOscillatorl(t,x,m0,r,w,m,f,k)
xdot=[x(2);-f*x(2)/(m0+m)-k*x(1)/(m0+m)+m0*r*cos(w*t)*w^2/(m0+m)]
错误:
??? Error using ==> plot
Vectors must be the same lengths.
Error in ==> LinearResponse at 22
plot(t,x3);
请各位大侠多多指教,谢谢!
[ 本帖最后由 ChaChing 于 2010-4-10 15:39 编辑 ] |