quadl求数值积分老出错

wangjizhe

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次，不方便用符号
望赐教

ChaChing

ref to 12F
常见的程序出错问题整理 (eight)
http://forum.vibunion.com/forum/thread-46001-1-1.html

wangjizhe

这是高人写的一篇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))
没有语法上的错误，但就是不对
所以我想我的也不是语法上的错误，但我不知道哪错了  注：不用符号积分

wangjizhe

我改用了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 编辑 ]

wangjizhe

我把z改为z(1),y改为y(1)，还是不对，版主帮我看看啥

我觉得quadl针对的是矢量，而我用的是标量，所以才出的错，但是我不知道怎样把它弄成矢量的形式，我也用符号积分了下，用了1分钟，但我得循环上千次呢，不可能叫我一直等下去吧

[ 本帖最后由 ChaChing 于 2009-5-17 18:53 编辑 ]

ChaChing

勿催帖, 现我在外头无matlab可试, 晚点有空试看看
但数值积分的东西, 我是来这里为回覆才学的, 以前没用过, 不见得对lz有帮助

wangjizhe

好的，真的非常感谢你

ChaChing

LZ的原始式子到底是什麽?
看了半天, 不明白:@)

[ 本帖最后由 ChaChing 于 2009-5-17 23:04 编辑 ]

wangjizhe

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)*[cos(B(k)) -sin(B(k)) 0 0]';
end
c=0.5/180*pi;
G=[1 0 1 0;0 1 0 1;0 0 1 0;0 0 0 1];
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)=[aa1,aa2,0,0]';
       A2(:,w)=[aa3,aa4,0,0]';
end

ChaChing

不知道是否因刚开完1~2hr的车, 看LZ的帖, 更晕:loveliness:
个人水平专业有限, 待高人路过

[ 本帖最后由 ChaChing 于 2009-5-17 23:35 编辑 ]

wangjizhe

明天详细给你说说，今天要熄灯了

wangjizhe

上面的代码不要看了，给你个简单的列子
function y=fun(t)
function f=f(x)
c1=[1,2,3,4];
c2=[1,0,x(1)-1,0;0,1,0,x(1)-1;0,0,1,0;0,0,0,1];
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.的错误，所以想找个方法解决这个问题。
帮我看看吧，版主，不然没法做了