帮忙看一下关于龙格库塔解微分方程的程序！（裂纹转子）

飞驰吧少年

我是刚入门做裂纹转子的，编程也是新手，从里到外都不太懂，麻烦各位大神帮忙看一下程序到底哪里出了问题？？
----------------------------------第一个m文件--------------------------------------------------------------
%jeffliewen
function dx=jeffliewen(t,x)
global w
m=32.1;
c=2100;
e=0.002;
k=2.5e7;
g=9.8;
s0=m*g/k;
omega=sqrt(k/m);
beta=0;
theta=pi;
k(1)=1-beta/pi+sin(4*beta)/(4*pi);%fa xiang gang du
k(2)=1-beta/pi+2*sin(2*beta)/(3*pi)-sin(4*beta)/(12*pi);%qie xiang gang du
k(3)=0;%ou he gang du
k(11)=k(1)*(cos(theta))^2+k(2)*(sin(theta))^2-k(3)*sin(2*theta);%x fang xiang gang du
k(22)=k(1)*(sin(theta))^2+k(2)*(cos(theta))^2+k(3)*sin(2*theta);%y fang xiang gang du
k(33)=(k(1)-k(2))*sin(theta)*cos(theta)+k(3)*cos(2*theta);%x,y fang xiang gang du
dx=[x(2);
    -x(2)*c/(m*(omega))-x(1)*k(11)/k-x(3)*k(33)/k+1+cos((theta)+(peta))*e*w^2/(s(0)*(omega)^2);
    x(4);
    -x(4)*c/(m*(omega))-x(3)*k(22)/k-x(1)*k(33)/k+sin((theta)+(peta))*e*w^2/(s(0)*(omega)^2)];

-------------------------------------------------第二个m文件-------------------------------------------------------------
%jeffliewentu
clear;
clc;
hold on;
global w
f=10:3:1200
for i=1:length(f)
    disp(f(i));
    w=f(i);
T=2*pi;
x(0)=[0.1;0;0.1;0];
tspan=0:T/100:100*T;
[t,x]=ode45('jeffliewen',tspan,x0);
x0=x(end,:);%ba yi ge zhou qi de chu zhi geng xin
plot(f(i),x(4000:100:end,5),'markersize',5);
xlabel ('转速w');
ylabel ('位移x1');
end


希望各位帮我指出一下错误，程序调试得在下要崩溃了。。。。
多谢~！！！

飞驰吧少年

我自己发现了第一处错：beta写成了peta……真是不细心。

yghit08

ode45，matlab不是有直接的命令么？
龙格库塔很简单。

飞驰吧少年

yghit08 发表于 2014-5-12 19:47
ode45，matlab不是有直接的命令么？
龙格库塔很简单。

对呀，我就是用的ode45，就是附上的那个程序，但是运行出来不行，老是有错误，我调试半天也不知道是哪里的问题，你能帮我看一下这个程序吗？谢谢~！

chybeyond

jeffliewen(t,x)，t好像在函数里没有用到，还有就是s(0)与s0,x0与x(0)

yghit08

好像是你的函数写错了，目前只能手机上网，不方便调试和查看

飞驰吧少年

yghit08 发表于 2014-5-12 21:01
好像是你的函数写错了，目前只能手机上网，不方便调试和查看

那可不可以麻烦你有空的时候上电脑帮我看看，这个对我满重要的，可是一直不能出图，我也不太明白到底是问题出在哪里。多谢~！

飞驰吧少年

chybeyond 发表于 2014-5-12 20:13
jeffliewen(t,x)，t好像在函数里没有用到，还有就是s(0)与s0,x0与x(0)

打不打括号是有区别的吗？那如果我改的话是统一改成打括号还是不打呢？？？

yghit08

你用的不是ode45这一类，是自己编程实现龙格库塔方法。建议查看帮助文档，另我两周多的时间里只能用手机

chybeyond

飞驰吧少年 发表于 2014-5-12 21:07
打不打括号是有区别的吗？那如果我改的话是统一改成打括号还是不打呢？？？

x(0)在matlab里边是错误的表达方式，按你的意思改成不打括号

飞驰吧少年

chybeyond 发表于 2014-5-12 21:38
x(0)在matlab里边是错误的表达方式，按你的意思改成不打括号

好的，十分感谢！！！但是如果不加括号的话在进行乘除的时候不会有歧义吗？

chybeyond

飞驰吧少年 发表于 2014-5-13 08:26
好的，十分感谢！！！但是如果不加括号的话在进行乘除的时候不会有歧义吗？

matlab变量命名规则：
1.变量名对大小写敏感
2.变量名第一个字符必须为英文字母，其长度不能超过31个字符；
3.变量名可以包含下连字符、数字，但不能包含空格符、标点。

飞驰吧少年

yghit08 发表于 2014-5-12 21:37
你用的不是ode45这一类，是自己编程实现龙格库塔方法。建议查看帮助文档，另我两周多的时间里只能用手机

非常感谢！！没关系，那以后就请教你一些小问题~
首先在ode45语句里调用函数时用@引出和用单引号括起来是一样的吗？
再有怎么查看帮助文档？
完全新手，多谢指教

飞驰吧少年

chybeyond 发表于 2014-5-13 08:42
matlab变量命名规则：
1.变量名对大小写敏感
2.变量名第一个字符必须为英文字母，其长度不能超过31个字 ...

谢谢！！！那我的k（11）什么之类的变量名也错了吧？我改改试试

chybeyond

飞驰吧少年 发表于 2014-5-13 10:36
谢谢！！！那我的k（11）什么之类的变量名也错了吧？我改改试试

k(11)这个没错，表示对向量k中的第十一个元素赋值，未赋值的元素默认为0，引用时直接使用k(11)即可，因为matlab中向量下标索引只能是正数或者逻辑数值，所以x(0)这种表达方式不对。