求助：一个求积分的程序

suntht

要绘制下面的函数随a的变化曲线，a从1到200取值。


程序如下：
fid=fopen('eta1.xls','w');
for i=linspace(1.0,200.0);
    a=i;

     q1=(a-x).*(x+1.0);
q2=sqrt(log(x+1.0));
   fun1=@(x)q1.*q2;
        q3=quadl(fun1,0,a);
      fun2=@(x)q1./q2;
     q4=quadl(fun2,0,a);  

eta=2.0*q3/q4;
    fprintf(fid,'%f %f\n',a,eta);
end

计算总是报错，请帮我看看问题出在哪儿？谢谢！

TBE_Legend

（1）matlab有个免费的package，叫 chebfun，很适合你的问题，超级强大。

（2）mmtc的结果供参考

winner245

你代码的错误主要在于匿名函数的定义部分。去掉文件读写，直接修改你的代码为（修改部分已经注释）：
A = linspace(1.0,200.0);
eta = zeros(size(A));
for n=1:length(A);
    a=A(n);
    q1=@(x)(a-x).*(x+1.0); % 函数定义 @(x)
    q2=@(x)sqrt(log(x+1.0));  % 函数定义 @(x)
    fun1=@(x)q1(x).*q2(x);  % 函数定义
    q3=quadl(fun1,0,a);
    fun2=@(x)q1(x)./q2(x); % 函数定义
    q4=quadl(fun2,0,a);
    eta(n)=2.0*q3/q4;
end
plot(A,eta)

另外，这类问题如果想高效计算的话，matlab提供了向量化积分解决方案（只要版本不低于2012a），效率为大幅度提升
a = linspace(1,200);
f1 = @(x) (a-x).*(x+1).*sqrt(log(x+1));
f2 = @(x) (a-x).*(x+1)./sqrt(log(x+1));
eta = 2*integral(@(x) f1(a*x).*a,0,1,'ArrayValued',true)./integral(@(x) f2(a*x).*a,0,1,'ArrayValued',true);
plot(a,eta)

suntht

TBE_Legend 发表于 2014-8-27 21:27
（1）matlab有个免费的package，叫 chebfun，很适合你的问题，超级强大。

（2）mmtc的结果供参考

谢谢推荐，chebfun很好用。
mmtc是什么？可以介绍一二吗？

suntht

winner245 发表于 2014-8-28 00:31
你代码的错误主要在于匿名函数的定义部分。去掉文件读写，直接修改你的代码为（修改部分已经注释）：
A = l ...

谢谢帮我改程序。因为我需要把几条曲线放到一张图上，想用画图软件处理数据，所以需要把结果输出文本文件，可以再麻烦你帮我把输出文件的代码补充上吗？谢谢！

winner245

suntht 发表于 2014-8-28 10:03
谢谢帮我改程序。因为我需要把几条曲线放到一张图上，想用画图软件处理数据，所以需要把结果输出文本文件 ...

我看你的代码是写 excel 文件，建议用 xlswrite 函数，查查帮助文件，照着例子写

suntht

winner245 发表于 2014-8-28 10:47
我看你的代码是写 excel 文件，建议用 xlswrite 函数，查查帮助文件，照着例子写 ...

写好了，谢谢帮忙。