quadl求数值积分老出错
function y=cv(x,p)global a;
global A1;
global F;
c=0.5*pi/180;
T=600;
GG=eye(4);
GG(1,3)=x(1)-T;
GG(2,4)=x(1)-T;
b=a(:,p);
y=2/c/c*b.'*GG*inv(F)*GG.'*A1(:,p);
*********************************************************************
a(4*600),A1(4*600),F(4*4)是在已求出的矩阵，我要对y求数值积分，但总出错
quadl('cv',600,60,1e-007,0,60)
??? Index exceeds matrix dimensions.
Error in ==> C:\MATLAB6p5\toolbox\matlab\funfun\quadl.m
On line 68==> if ~isfinite(y(13))，
因为我要积分600次，不方便用符号
望赐教
回复 楼主 wangjizhe 的帖子
ref to 12F常见的程序出错问题整理 (eight)
http://forum.vibunion.com/forum/thread-46001-1-1.html
回复 沙发 ChaChing 的帖子
这是高人写的一篇matlab编程小技巧http://forum.vibunion.com/thread-2588-1-1.html
a = rand;
b = rand + 1;
K = rand;
y = subs('sin(K*x)/x', 'K', sym(K, 'd') )
quad( inline(vectorize(char(y))), a, b);
第二，注意vectorize函数。sym对象重载的char函数会把数组的乘除(.*, ./, .^)化成矩阵的乘除(*, /, ^)，vectorize会将字符串里的*,/,^全化成.*, ./, .^，如果不这样做，quad函数将会出错，这个问题曾经困扰了我很长的时间：
quad( inline(char(y)), a, b)
??? Index exceeds matrix dimensions
Error in ==> D:\MATLAB6p5\toolbox\matlab\funfun\quad.m
On line 67 ==> if ~isfinite(y(7))
没有语法上的错误，但就是不对
所以我想我的也不是语法上的错误，但我不知道哪错了注：不用符号积分 我改用了7.1的，rocwoods 的列子对了，可我的还是算不出来，不知道为什么呀，麻烦帮我看看好吗
function y=cv(z)
global a A1 F;
T=600;
function f=f(x)
c=0.5*pi/180;
T=600;
GG=eye(4);
GG(1,3)=x(1)-T;
GG(2,4)=x(1)-T;
b=a(:,z);
f=2/c/c*b.'*GG*inv(F)*GG.'*A1(:,z);
end
y=quadl(@f,T,z);
end
cv(20)
??? Attempted to access y(13); index out of bounds because numel(y)=1.
Error in ==> quadl at 72
if ~isfinite(y(13))
Error in ==> cv at 13
y=quadl(@f,T,z);
[ 本帖最后由 ChaChing 于 2009-5-17 19:01 编辑 ] 我把z改为z(1),y改为y(1)，还是不对，版主帮我看看啥
我觉得quadl针对的是矢量，而我用的是标量，所以才出的错，但是我不知道怎样把它弄成矢量的形式，我也用符号积分了下，用了1分钟，但我得循环上千次呢，不可能叫我一直等下去吧
[ 本帖最后由 ChaChing 于 2009-5-17 18:53 编辑 ] 勿催帖, 现我在外头无matlab可试, 晚点有空试看看
但数值积分的东西, 我是来这里为回覆才学的, 以前没用过, 不见得对lz有帮助
回复 6楼 ChaChing 的帖子
好的，真的非常感谢你 LZ的原始式子到底是什麽?看了半天, 不明白:@)
[ 本帖最后由 ChaChing 于 2009-5-17 23:04 编辑 ]
回复 8楼 ChaChing 的帖子
a(4*600),A1(4*600),F(4*4)是在已求出的矩阵T=600;
U=eye(1,T);
Xo=zeros(1,T);
Yo=zeros(1,T);
Vo=3.0;
global a;
global A1;
global F;
for k=1:T-1
Xo(k+1)=Xo(k)+Vo*sin(U(k));
Yo(k+1)=Yo(k)+Vo*cos(U(k));
end
Xt=zeros(1,T);
Yt=1e+4*eye(1,T);
Vt=5.0;
for k=1:T-1
Xt(k+1)=Xt(k)+Vt*cos(pi/3);
Yt(k+1)=Yt(k)-Vt*sin(pi/3);
end
rx=zeros(1,T);
ry=zeros(1,T);
r=zeros(1,T);
for k=1:T
rx(k)=Xt(k)-Xo(k);
ry(k)=Yt(k)-Yo(k);
r(k)=sqrt(rx(k)^2+ry(k)^2);
end
B=atan2(rx,ry);
a=zeros(4,T);
for k=1:T
a(:,k)=1/r(k)*';
end
c=0.5/180*pi;
G=;
F=zeros(4);
for k=1:T
F=G.'*F*G+1/c/c*a(:,k)*a(:,k).';
end
for w=1:T
Xtt=Xt(w);
Xoo=Xo(w);
Ytt=Yt(w);
Yoo=Yo(w);
aa1 =-1/2/(Xtt^2-2*Xtt*Xoo+Xoo^2+Ytt^2-2*Ytt*Yoo+Yoo^2)^(3/2)/(1+(Xtt-Xoo)^2/(Ytt-Yoo)^2)^(1/2)*(-2*Xtt+2*Xoo)+1/(Xtt^2-2*Xtt*Xoo+Xoo^2+Ytt^2-2*Ytt*Yoo+Yoo^2)^(1/2)/(1+(Xtt-Xoo)^2/(Ytt-Yoo)^2)^(3/2)*(Xtt-Xoo)/(Ytt-Yoo)^2;
aa2 =1/2/(Xtt^2-2*Xtt*Xoo+Xoo^2+Ytt^2-2*Ytt*Yoo+Yoo^2)^(3/2)*(Xtt-Xoo)/(Ytt-Yoo)/(1+(Xtt-Xoo)^2/(Ytt-Yoo)^2)^(1/2)*(-2*Xtt+2*Xoo)+1/(Xtt^2-2*Xtt*Xoo+Xoo^2+Ytt^2-2*Ytt*Yoo+Yoo^2)^(1/2)/(Ytt-Yoo)/(1+(Xtt-Xoo)^2/(Ytt-Yoo)^2)^(1/2)-1/(Xtt^2-2*Xtt*Xoo+Xoo^2+Ytt^2-2*Ytt*Yoo+Yoo^2)^(1/2)*(Xtt-Xoo)^2/(Ytt-Yoo)^3/(1+(Xtt-Xoo)^2/(Ytt-Yoo)^2)^(3/2);
aa3 =-1/2/(Xtt^2-2*Xtt*Xoo+Xoo^2+Ytt^2-2*Ytt*Yoo+Yoo^2)^(3/2)/(1+(Xtt-Xoo)^2/(Ytt-Yoo)^2)^(1/2)*(-2*Ytt+2*Yoo)-1/(Xtt^2-2*Xtt*Xoo+Xoo^2+Ytt^2-2*Ytt*Yoo+Yoo^2)^(1/2)*(Xtt-Xoo)^2/(Ytt-Yoo)^3/(1+(Xtt-Xoo)^2/(Ytt-Yoo)^2)^(3/2);
aa4 =1/2/(Xtt^2-2*Xtt*Xoo+Xoo^2+Ytt^2-2*Ytt*Yoo+Yoo^2)^(3/2)*(Xtt-Xoo)/(Ytt-Yoo)/(1+(Xtt-Xoo)^2/(Ytt-Yoo)^2)^(1/2)*(-2*Ytt+2*Yoo)-1/(Xtt^2-2*Xtt*Xoo+Xoo^2+Ytt^2-2*Ytt*Yoo+Yoo^2)^(1/2)*(Xtt-Xoo)/(Ytt-Yoo)^2/(1+(Xtt-Xoo)^2/(Ytt-Yoo)^2)^(1/2)+1/(Xtt^2-2*Xtt*Xoo+Xoo^2+Ytt^2-2*Ytt*Yoo+Yoo^2)^(1/2)*(Xtt-Xoo)^3/(Ytt-Yoo)^4/(1+(Xtt-Xoo)^2/(Ytt-Yoo)^2)^(3/2);
A1(:,w)=';
A2(:,w)=';
end 不知道是否因刚开完1~2hr的车, 看LZ的帖, 更晕:loveliness:
个人水平专业有限, 待高人路过
[ 本帖最后由 ChaChing 于 2009-5-17 23:35 编辑 ]
回复 6楼 ChaChing 的帖子
明天详细给你说说，今天要熄灯了回复 10楼 ChaChing 的帖子
上面的代码不要看了，给你个简单的列子function y=fun(t)
function f=f(x)
c1=;
c2=;
f=c1*c2*c1';
end
y=quadl(@f,0,t);
end
***************************
问题的症结：quadl是变步长辛普森求积，y=quadl(@f,0,t);中的变量是x,通过调式，x 有13个值(因为是变步长嘛)，而c2中我用的是x(1),所以出错；简单点所就是quadl针对的是矢量，而我把它标量化了.但是因为c2中有个x(1)-1,如果我改为x-1;会出现??? Error using ==> vertcat
All rows in the bracketed expression must have the same number of columns.的错误，所以想找个方法解决这个问题。
帮我看看吧，版主，不然没法做了
页:
[1]