声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1444|回复: 1

[编程技巧] 一个常微分的方程组调用龙格库塔四阶解

[复制链接]
发表于 2011-7-22 12:38 | 显示全部楼层 |阅读模式

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

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

x
Function dx=dpeq_one(t,x)
dx=zeros(3,1);
dx(1)=(-3.43771603048094e-05*x(2)^6-0.000132921305442425*x(2)^5+0.00415935962144218*x(2)^4+0.0233126081359259*x(2)^3-...
0.0732414070985704*x(2)^2-1.05403016904226*x(2)-5.85820689312774)*x(3)+(-4.35789199696095e-05*x(2)^6-2.69424808482508e-05*x(2)^5+...
0.00635716041316867*x(2)^4+0.0185421269112098*x(2)^3-0.220192402342151*x(2)^2-1.38426845614949*x(2)-6.76278220704981)*x(1)+...
(1.13790540618933e-06*x(2)^6+1.41177469687368e-05*x(2)^5-5.94593881351547e-05*x(2)^4-0.00150168194443906*x(2)^3-...
0.00538960610161732*x(2)^2+0.0234979236508000*x(2)+0.172445454319882)*x(3)*x(1)+...
(7.15581608696399e-06*x(2)^6+4.42405268444341e-06*x(2)^5-0.00104386870495382*x(2)^4-0.00304468422187366*x(2)^3+...
0.0361563879051154*x(2)^2+0.227301881141133*x(2)+1.11047326880949)*sin(108*pi*t);%+...(7.98181866015145e-08*x(2)^6+4.63704589950561e-07*x(2)^5-8.61063595839847e-06*x(2)^4-6.81446324712728e-05*x(2)^3+9.82619193232414e-05*x(2)^2+...0.00132096154647624*x(2)+0.0193684391761411)*x(1)*x(2);
dx(2)=x3;
dx(3)=((-2.34890000000000e-05*x(2)^4+4.33260000000000e-05*x(2)^3-0.00328580000000000*x(2)^2+0.00276810000000000*x(2))*x(2)-1.283*x(3)+...
    (-1.01490000000000e-06*x(2)^6-3.60750000000000e-06*x(2)^5+0.000381580000000000*x(2)^4+0.000832520000000000*x(2)^3-0.0616820000000000*x(2)^2-...
    0.0575550000000000*x(2)+5.21430000000000)*x(1)+(7.51529691386400e-21*x(2)^6-5.05625999999975e-07*x(2)^5-1.56135000000007e-05*x(2)^4+...
    7.85119999999977e-05*x(2)^3+0.00277059000000001*x(2)^2+0.00208740000000004*x(2)-0.0742719999999999)*x(1)^2)/0.015593;
x=[x(1),x(2),x(3)];
使用下来函数来解:
clear
clc
%获得时间序列
[t,x]=ode45(@dpeq_one,[0,400],[0;0;0]);
subplot(2,1,1)
plot(t,x(:,1))
subplot(2,1,2)
plot(t,x(:,3))

h2=figure(2);
subplot(2,1,1)
plot(t,x(:,2))

h3=figure(3);
plot(x(:,1),x(:,3))

h4=figure(4);
plot(t,x(:,1),'r',t,x(:,3),'b')  %或者用plot(t,y(:,1),'-',t,y(:,2),'-.')
出现以下错误:
Warning: feval on script names will not work, or may work differently,
in a future version of MATLAB.  To make your code insensitive to any change
and to suppress this warning message:
- Either change the script to a function.
- Or use eval instead of feval.
The script file in question is dpeq_one.
> In funfun\private\odearguments at 110
  In ode45 at 173
??? Error using ==> feval
Attempt to execute SCRIPT dpeq_one as a function:
C:\Users\asus\Documents\MATLAB\dpeq_one.m

Error in ==> odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
回复
分享到:

使用道具 举报

发表于 2011-7-22 16:01 | 显示全部楼层
...Attempt to execute SCRIPT dpeq_one as a function:...

由此判断函数定义有问题!
matlab是有大小区分的, 试试将Function改成function
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-5-6 17:01 , Processed in 0.056700 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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