声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 5153|回复: 2

[编程技巧] 请帮忙修改一下四阶龙格库塔法解微分方程的matlab程序

[复制链接]
发表于 2007-5-16 13:35 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
请帮忙编个四阶龙格库塔法解非线性微分方程的程序给我看看,用matlab语言的,方程是
          xdot(2)=(-(1/tau1+K*tau2*(cosx(1))/tau1)*x(2)-K*(sinx(1))/tau1+Deltaomegao/tau1);
         xdot(1)=x(2);
其中xdot(1)就是x1求一阶导,xdot(2)就是x2求一阶导,Deltaomegao,tau1,tau2,K是要给定初值的常量
我把其中的cosx(1),sinx(1)用泰勒公式展开了前三项然后用ode45编了程序,但是老师说这样不精确,也不专业,要我用龙格库塔法,不会,所以请大家帮帮忙,看看怎么帮我改改我的程序。
本人的程序如下:
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=[1,2,3,4];
tspan=linspace(0,15,1500);
options=odeset('Abstol',[1e-6,1e-6]);
for i=1:4
     [t,x]=ode45('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:'(

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

 楼主| 发表于 2007-5-16 14:28 | 显示全部楼层

请教一个用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=[1,2,3,4,5,6,7,8,9,10];
tspan=linspace(0,15,1500);
options=odeset('Abstol',[1e-6,1e-6]);
for i=1:10
     [t,x]=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:39 | 显示全部楼层

请帮看一下程序运行中的错误怎么改?

将下程序保存为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);
程序ode4.m  四阶龙格库塔法
function Y = ode4(odefun,tspan,y0,varargin)
if ~isnumeric(tspan)
error('TSPAN should be a vector of integration steps.');
end
if ~isnumeric(y0)
error('Y0 should be a vector of initial conditions.');
end
h = diff(tspan);
if any(sign(h(1))*h <= 0)
error('Entries of TSPAN are not in order.')
end
try
f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
error(msg);
end
y0 = y0(:); % Make a column vector.
if ~isequal(size(y0),size(f0))
error('Inconsistent sizes of Y0 and f(t0,y0).');
end
neq = length(y0);
N = length(tspan);
Y = zeros(neq,N);
F = zeros(neq,4);
Y(:,1) = y0;
for i = 2:N
ti = tspan(i-1);
hi = h(i-1);
yi = Y(:,i-1);
F(:,1) = feval(odefun,ti,yi,varargin{:});
F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});
F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2),varargin{:});
F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3),varargin{:});
Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));
end
Y = Y.';
主程序:
tau1=5;tau2=1;K=12;Deltaomegao=[1,2,3,4,5,6,7,8,9,10];
tspan=linspace(0,15,1500);
options=odeset('Abstol',[1e-6,1e-6]);
for i=1:10
     [t,x]=ode4('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
本人找到一个四阶龙格库塔法的程序,但是在主程序中调用时出现如下错误,:@L 请教一下怎么修改?:@L
??? Error using ==> ode4
Too many output arguments.
Error in ==> TwoPLLshzhjiyituduxi at 5
     [t,x]=ode4('TwoPLL2',tspan,[-pi,Deltaomegao(i)]',options,tau1,tau2,K,Deltaomegao(i));
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-24 23:24 , Processed in 0.065172 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表