Question

求助，用matlab积分求面积问题，快要上交，求速度~~~~

被积函数是对数概率密度函数，用simpson法数值计算出[0,2000]内的面积，这是我写的，第一个函数是simpson法的定义，第二个函数是对数正态分布的概率密度函数，但是运行出错，我看不出那里出错了，求助。。。。。。速度。。。。。
function s=MySimp(x,h)
n=length(x);
ys=f(x(1));
yf=f(x(n));
ymid=0;
yin=0;
for i=1:n-1
ymid=ymid+4*f((x(i)+x(i+1))/2);
end
for i=2:n-1
yin=yin+2*f(x(i));
end
s=(h/6)*(ys+ymid+yin+yf);

function y=f(x)
y=(1/x*0.891*sqrt(2*pi))*exp((-(logx-7.59)^2)/2*0.7938);
a=0;
b=2000;
n=8;
h=(b-a)/8;
x=[0:h:2000];
ssimp=MySimp(x,h);
ssimp

Answer 1

S=（L2的曲线积分－L1的曲线积分）
用三次样条插值和复化辛普森公式计算：
例如
x0=[0.9 1.3 1.7 2.1 2.6 3 3.2 3.3 3.5 4.0 4.6 5.0 5.5 6.0 6.3 6.6 7.0 7.3 7.8 8.1 8.4 8.8 9.3 9.6 10.0 10.2 10.6 10.8 11.1 11.5 11.7 12.0 12.5 12.9 13.6 14.0 14.5 15.1 15.6 16.1 16.7 17.0 17.5 17.6 17.8];

y1=[6.5 6.0 5.6 5.6 5.0 4.5 4.2 4.1 4.1 3.9 3.6 3.0 3.4 3.5 3.5 3.5 3.4 3.1 3.2 3.4 3.3 3.3 3.8 3.2 2.8 2.5 2.2 2.4 2.6 2.0 2.6 2.9 2.9 3.4 1.5 1.1 1.0 0.7 1.9 1.8 2.2 2.6 2.6 2.7 3.0]; %已知节点（x0,y1）

y2=[6.5 7.3 7.2 7.3 7.7 7.3 7.6 9.9 10.1 10.5 10.2 9.8 9.7 8.8 8.6 8.2 7.7 7.5 7.7 8.2 8.9 8.2 9.7 9.7 9.7 10.0 10.1 10.4 10.3 10.3 10.5 10.5 10.4 10.2 9.6 9.7 8.9 7.9 7.6 5.9 5.8 3.4 3.3 3.2 3.0]; %已知节点（x0,y2）

h=0.1; %产生插值点的步长

x=0.9:h:17.8; %产生插值点x

y01=interp1(x0,y1,x); %计算L1的分段线性插值

y02=interp1(x0,y2,x); %计算L2的分段线性插值

plot(x,y01,'k',x,y02,'r') %分段线性插值作图

z1=trapz(y01)*h %用梯形公式计算L1 的积分
>> y021=[y02(2:2:k-1)];

>> s021=sum(y021);

>> y022=[y02(3:2:k-1)];

>> s022=sum(y022);

>> z2=(y02(1)+y02(k)+4*s021+2*s022)*h/3 %用辛普森公式计算L2的积分
S=(z2-z1) %代入面积推导公式

>> x0=0:1:24;
y0=[0 0 1.01 1.18 0 1.39 1.97 1.33 2.08 1.22 1.09 0 1.01 0 0.9 1.18 1.13 1.05 1.58 0.94 0 0 1.18 1.01 0.9];
h=0.1; %产生插值点的步长
x=0:h:24; %产生插值点x
y01=interp1(x0,y0,x); %计算L1的分段线性插值
plot(x,y01,'k') %分段线性插值作图
z1=trapz(y01)*h %用梯形公式计算L1 的积分

z1 =

21.7000

样条插值
x0=0:1:24;
y0=[0 0 1.01 1.18 0 1.39 1.97 1.33 2.08 1.22 1.09 0 1.01 0 0.9 1.18 1.13 1.05 1.58 0.94 0 0 1.18 1.01 0.9];
h=0.1; %产生插值点的步长
x=0:h:24; %产生插值点x
y01=spline(x0,y0,x); %计算L1样条插值
plot(x,y01,'k') %分段线性插值作图
z1=trapz(y01)*h %用梯形公式计算L1 的积分

z1 =

21.5442

Answer 2

function s=MySimp(a,b,k)
h=(b-a)/k;
x=[a:h:b];
n=length(x);
ys=f(x(2));
yf=f(x(n));
ymid=0;
yin=0;
for i=1:n-1
ymid=ymid+4*f((x(i)+x(i+1))/2);
end
for i=2:n-1
yin=yin+2*f(x(i));
end
s=(h/6)*(ys+ymid+yin+yf);

function y=f(x)
y=(1/x*0.891*sqrt(2*pi))*exp((-(log(x)-7.59)^2)/2*0.7938);
a=0;
n=8;
b=2000;
ssimp=MySimp(a,b,n)

ssimp =

3.2412

追问

为什么运行之后还是有错，运行后出现这样的结果。。。。望多指教一下。。。。谢谢。。。。
ssimp=MySimp(a,b,n)
??? Error using ==> MySimp
Too many input arguments.

追答

你的函数改了吗，我运行没问题

Answer 3

第二个函数第一行：
logx语法错误吧