Question

求MATLAB大神帮我看看程序哪有问题，谢谢！

搞了一下午了，拜托大家帮忙！
function g=f()
clear all;
t=0.1;
R=0.5;
a=2; %D的取值区间上限;
b=0.04; %D的步长;
c1=0; %O取值区间下限;
c2=10; %O取值区间上限;
d=0.1; %O的步长;
P21=1;r21=0.1;
syms y1;
g=zeros(((c1-c2)/d+1),(a/b*2+1));
n=0;
for O=c1:d:c2
if (O-fix(O)==0)
clear maplemex;
end
n=n+1;
m=0;
for x=-a:b:a
m=m+1;
Rs=(O^2*(R+2*r21))/(4*(x+0.002*y1)^2+(R+2*r21)^2);
D=4*Rs^2+R^2+4*Rs*P21+2*P21*R+P21^2;
s1=0.5*(-(2*Rs+P21+R)+sqrt(D));
s2=0.5*(-(2*Rs+P21+R)-sqrt(D));
r11=Rs*(Rs+P21)/sqrt(D)*(1/(Rs+s2)*(exp(-Rs*t)-exp(s2*t))-1/(Rs+s1)*(exp(-Rs*t)-exp(s1*t)))+exp(-Rs*t); %ρ11
r22=Rs/sqrt(D)*(exp(s1*t)-exp(s2*t));
r12=0.5*O*(R/2+r21)*(r11-r22)/((x+0.002*y1)^2+(R/2+r21)^2);
f=2.4*10^(-8)*y1.^2.*exp(-4.83*10^(-6)*y1^2)*r12; %得到被积函数
f=@(y)subs(f,y1,y);
g(n,m)=quad(f,0,100);
end
end
y=c1:d:c2;
x=-a:b:a;
Z=g;
surf(x,y,Z);
end

这个程序运行完屏幕上显示
>> f;
??? Attempted to access y(7); index out of bounds because numel(y)=1.

Error in ==> quad at 71
if ~isfinite(y(7))

Error in ==> f at 31
g(n,m)=quad(f,0,100);
每循环一次，y就会被赋值一次，怎么会出现y长度不够的问题呢？

Answer 1

我把步长改了改, 因为不然电脑运行太慢不出结果, 所以改长了, 但是没有其他影响. 下面的代码是没改的.

t=0.1;
R=0.5;
a=2; %D的取值区间上限;
b=0.04; %D的步长, 我取的0.1; 
c1=0; %O取值区间下限;
c2=10; %O取值区间上限;
d=0.1; %O的步长, 我取的0.5;
P21=1;r21=0.1;
syms y1;
g=zeros(((c1-c2)/d+1),(a/b*2+1));
n=0;
for O=c1:d:c2
if (O-fix(O)==0)
clear maplemex;
end
n=n+1;
m=0;
for x=-a:b:a
m=m+1;
Rs=(O^2*(R+2*r21))/(4*(x+0.002*y1)^2+(R+2*r21)^2); 
D=4*Rs^2+R^2+4*Rs*P21+2*P21*R+P21^2; 
s1=0.5*(-(2*Rs+P21+R)+sqrt(D)); 
s2=0.5*(-(2*Rs+P21+R)-sqrt(D));
r11=Rs*(Rs+P21)/sqrt(D)*(1/(Rs+s2)*(exp(-Rs*t)-exp(s2*t))-1/(Rs+s1)*(exp(-Rs*t)-exp(s1*t)))+exp(-Rs*t); %ρ11
r22=Rs/sqrt(D)*(exp(s1*t)-exp(s2*t)); 
r12=0.5*O*(R/2+r21)*(r11-r22)/((x+0.002*y1)^2+(R/2+r21)^2); 
f=2.4*10^(-8)*y1.^2.*exp(-4.83*10^(-6).*y1.^2).*r12; %得到被积函数
f=@(y1)subs(f)+0.*y1;  % 这里加了一项0, 就能运行了, 结果肯定不受影响.
g(n,m)=quad(f,0,100); 
end
end
y=c1:d:c2;
x=-a:b:a;
Z=g;
surf(x,y,Z);

图像是精简之后的, 不过你画出来应该是类似的, 仅仅是更精细:

追问

确实可以啦！为什么必须要加“0.*y1”这一项呢？改之前的代码到底哪里有问题呢？

追答

那个quad函数的确有问题, 包括quadl函数也一样, 被积函数只要是常数就有这问题, 你可以试试. 不知道为什么. 加上一项0就行了.

追问

哦，那就是这个quad函数原本就有的bug了？万分感谢大神指点！

追答

应该是, 因为常值函数不含自变量, 可能就出错了, 我猜是这样吧. quad函数可以调试, 但是我没分析里边的代码. 反正能出来结果, 不求甚解了.

Answer 2

原因是f是你的函数的名字，又是你函数里变量的名字，导致冲突。你把你的函数改个名字就可以了。

追问

不是这里的问题哦，已经解决啦，多谢！

追答

你问题里的错误就是我说的原因导致的，另外quad没有问题，我也没有加什么0.*，只是把你程序里的f换成了ff就可以出来图，是这样的。

回你只是因为觉得你这样写程序用和函数名相同的名字取变量名会给自己造成麻烦，以后要注意。