insects 发表于 2007-5-16 14:31

请教一个用matlab解微分方程的小问题?

用matlab解微分方程组的时候,方程组中有sin,cos函数怎么办?怎么编呢?要不要先用泰勒展开先呢,ode45能不能直接解这样的方程组啊?
例:
function xdot=0000001(t,x,dummy,tau1,tau2,K,Deltaomegao)
         xdot=zeros(2,1);
         xdot(2)=(-(1/tau1+K*tau2*sin(x(1))/tau1)*x(2)-K*cos(x(1))/tau1+Deltaomegao/tau1);
         xdot(1)=x(2);
还是要展开先呢?将下程序保存为TwoPLL2.m
function xdot=TwoPLL2(t,x,dummy,tau1,tau2,K,Deltaomegao)
         xdot=zeros(2,1);
         xdot(2)=(-(1/tau1+K*tau2*(1-x(1)^2/2+x(1)^4/24)/tau1)*x(2)-K*(x(1)-x(1)^3/6+x(1)^5/120)/tau1+Deltaomegao/tau1);
         xdot(1)=x(2);
不过自己展开的话精确度就小了很多,因为自己展开不有写出很多项,只能展开前几项啊,能解决这个问题吗?
主程序:
tau1=5;tau2=1;K=12;Deltaomegao=;
tspan=linspace(0,15,1500);
options=odeset('Abstol',);
for i=1:10
   =ode15s('TwoPLL2',tspan,[-pi,Deltaomegao(i)]',options,tau1,tau2,K,Deltaomegao(i));
   figure(1);plot(t,x(:,1));       hold on; title('Phase difference as a function of Time');ylabel('θe(t)/rad');xlabel('t/s');
   figure(2);plot(t,x(:,2));       hold on; title('Frequency difference as a function of Time');ylabel('(dθe(t)/dt)/(rad?sˉ1)');xlabel('t/s');
   figure(3);plot(x(:,1),x(:,2));hold on; title('Phase portrait');ylabel('(dθe(t)/dt)/(rad?sˉ1)');xlabel('/rad');
end
可以运行一下,没有错误的,我有运行过了的。:handshake

咕噜噜 发表于 2007-5-16 14:37

回复 #1 insects 的帖子

能不展开就好还是不要展开的好,ode45对于通常的含有正弦余弦的函数还是可以求解的

gghhjj 发表于 2007-5-17 07:58

无需展开
页: [1]
查看完整版本: 请教一个用matlab解微分方程的小问题?