请教一个用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
回复 #1 insects 的帖子
能不展开就好还是不要展开的好,ode45对于通常的含有正弦余弦的函数还是可以求解的 无需展开
页:
[1]