求助关于求导的问题

蜜游儿

我用ode45求解了一个二阶线性非齐次方程，得到x1和x2=dx1/dt与时间t 的关系图，可是我想还得到dx2/dt与时间t 的关系图该怎么办呢？
ode方程是这样的：dx1/dt=x2；dx2/dt=f*x2-k*x1(f,k是已知量)

houxiaoyang

你用ode45返回的应该有X1与X2的值，再加一列X3，用其表示dx2/dt

ChaChing

除LS说的方式, 试试diff是否可行!?

蜜游儿

谢谢你啦！

蜜游儿

谢谢你啦！

蜜游儿

我用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 编辑 ]

houxiaoyang

你不是返回个X么，里面有X1,X2，你再加一列X3，X3可以用X2,X1表示吧，这样就可以了，我以前也遇到过和LZ一样的问题，就是这样解决的，虽然有点麻烦

maigicku

从错误看t与x3的长度不相同。。检查一下t和x3

蜜游儿

是的 这个问题已经解决了，就是这个方法，谢谢