带2参数的数值积分

yxlnbu

大家好，我有一个函数'2*(1600-40*r*cos(x))/(r^2-80*r*cos(x)+1600+z^2)^1.5，其中r和z是两个参数，r=【40：45】
z=【-4.5：4.5】我要对x进行从0到pi的数值积分，可是从workspace中看到每次积分结果都是一样，想请教各位大师了。在此先谢过了
clear all
clc
rr=linspace(40.00,45.00,100);%产生r=40到45的数组
zz=linspace(-4.5,4.5,6);%产生每匝的z位移从-4.5到4.5
b=zeros(100,8);
for j=1:100             %首先赋值使公式仅为x的方程，利用Simpson积分将在r的每一处b（r）求出，然后拟合为5次的多项式
    r=rr(j);            %对r赋值，使积分式中r为数值
    for i=1:6           
        z=zz(i);        %对z赋值，使积分式z1为数值
        f=strcat('2*(1600-40*',num2str(r),'*cos(x))/(',num2str(r),'^2-80*',num2str(r),'*cos(x)+1600+',num2str(z),'^2)^1.5');
        b(j,i)=quad(inline('f'),0,pi);%simpos数值积分，将产生的值放入到b（j，i）中
        b(j,7)=b(j,7)+b(j,i);%将每匝在试样环处的系数累计至第7列中
        b(j,8)=r;               %将r值放入第8列中，方便拟合
    end
end

zhouyang664

function  b = odeFun( )

close all ;
clear all ;
clc ;
format long ;

rr = linspace(40.00 , 45.00 , 100) ;          %产生r=40到45的数组
zz = linspace(-4.5 , 4.5 , 6) ;               %产生每匝的z位移从-4.5到4.5
b = zeros(100 , 8) ;

for j = 1 : 100                               %首先赋值使公式仅为x的方程，利用Simpson积分将在r的每一处b（r）求出，然后拟合为5次的多项式
    r = rr (j) ;                              %对r赋值，使积分式中r为数值
    for i = 1 : 6           
        z = zz(i) ;                           %对z赋值，使积分式z1为数值
        b(j , i) =  quad(@(x) 2*(1600 -40*r*cos(x))./(r^2-80*40*cos(x)+1600+z^2).^1.5 , 0 , pi) ;
                                              %simpos数值积分，将产生的值放入到b（j，i）中
        b(j , 7) = b(j , 7) + b(j , i) ;      %将每匝在试样环处的系数累计至第7列中
        b(j , 8) = r ;                        %将r值放入第8列中，方便拟合
    end
end
end

输出b为你要的数据 ！

yxlnbu

回复 2 # zhouyang664 的帖子

谢谢您，在这个问题上花了好长时间不得其解，谢谢您的帮助。能解释下为什么inline函数不行么？有时候就说输入不够之类的。

ChaChing

回复 3 # yxlnbu 的帖子

LZ比較下兩者差異, 應該即可體會
f='2*x+3';
ff=inline('f')
ff=inline(f)

yxlnbu

回复 4 # ChaChing 的帖子

谢谢，刚才我试了下，一个是以x为变量，另一个是以f为变量的，请问是这样理解么？另外，我把上面的程序去掉引号后，程序计算就出现了奇异，而应用zhouyang664的程序的时候就没有出现这个问题。但把rr和zz值缩小十倍之后也出现了奇异
Warning: Maximum function count exceeded; singularity likely.

happy

yxlnbu 发表于 2011-3-11 10:32

                               
登录/注册后可看大图

回复 4 # ChaChing 的帖子

谢谢，刚才我试了下，一个是以x为变量，另一个是以f为变量的，请问是这样理解么 ...

没看懂你的问题
直接把你修改后而出现问题的代码贴上来，是查找问题的最直接方法

yxlnbu

回复 6 # happy 的帖子
你好，问题是这样的，我把inline中的引号去掉（红色部分），尝试一下，然后就出现了奇异
出现的警告是这样的（总共出现上百次）：Warning: Maximum function count exceeded; singularity likely.
> In quad at 110
  In Untitled at 11

clear all
clc
rr=linspace(40.00,45.00,100);%产生r=40到45的数组
zz=linspace(-4.5,4.5,6);%产生每匝的z位移从-4.5到4.5
b=zeros(100,8);
for j=1:100             %首先赋值使公式仅为x的方程，利用Simpson积分将在r的每一处b（r）求出，然后拟合为5次的多项式
    r=rr(j);            %对r赋值，使积分式中r为数值
    for i=1:6           
        z=zz(i);        %对z赋值，使积分式z1为数值
        f=strcat('2.*(1600-40.*',num2str(r),'.*cos(x))./(',num2str(r),'.^2-80.*',num2str(r),'.*cos(x)+1600+',num2str(z),'.^2).^1.5');
        b(j,i)=quad(inline(f),0,pi);%simpos数值积分，将产生的值放入到b（j，i）中,红色是改掉的
        b(j,7)=b(j,7)+b(j,i);%将每匝在试样环处的系数累计至第7列中
        b(j,8)=r;               %将r值放入第8列中，方便拟合
    end
end

zhouyang664

楼主可以自己用断点调试查一下问题的所在啊！

TNC

q = quad(fun,a,b) tries toapproximate the integral of function fun from a to b towithin an error of 1e-6 using recursive adaptiveSimpson quadrature. fun is a function handle.

ChaChing

回复 7 # yxlnbu 的帖子

為何兩者一個會出警訊, 另一個則不會?

Warning: Maximum function count exceeded; singularity likely.

兩者的差異僅是函數的使用方式不同而已, 頂多速度/效率不同, 不應該有如此差異吧!?

f=strcat('2*(1600-40*',num2str(r),'*cos(x))./(',num2str(r),'^2-80*',num2str(r),'*cos(x)+1600+(',num2str(z),')^2).^1.5');
b(j,i)=quad(inline(f),0,pi);

b(j,i) =  quad(@(x) 2*(1600 -40*r*cos(x))./(r^2-80*40*cos(x)+1600+z^2).^1.5 , 0 , pi) ;

ChaChing

软件的警讯是不会错的! 也就是说函数有奇异, 但会因使用inline或@而不同吗?
我看楼主也是在休周末了, 就留下问题, 给有兴趣自我提升的人试试吧!

yxlnbu

回复 11 # ChaChing 的帖子

谢谢各位，我是回家了几天，没有上论坛了。今天刚回来。
函数中有改变的就是@和inline函数调用，还有的就是‘num2str’ 函数的调用。大家可以试试。

ChaChing

都忘了这了! 我想没人想试看看赚威望了
看看差异!!

f=strcat('2*(1600-40*',num2str(r),'*cos(x))./(',num2str(r),'^2-80*',num2str(r),'*cos(x)+1600+',num2str(z),'^2).^1.5');

f=strcat('2*(1600-40*',num2str(r),'*cos(x))./(',num2str(r),'^2-80*',num2str(r),'*cos(x)+1600+(',num2str(z),')^2).^1.5');

qibbxxt

回复 13 # ChaChing 的帖子

按说Matlab不会出现这种错误的，其实这两种写法，函数是不一样的

1. f=strcat('2*(1600-40*',num2str(r),'*cos(x))/(',num2str(r),'^2-80*',num2str(r),'*cos(x)+1600+',num2str(z),'^2)^1.5');

复制代码

这种写法，当z是负值的时候，就是

1. 2*(1600-40*40*cos(x))/(40^2-80*40*cos(x)+1600+-4.5^2)^1.5

复制代码

而后面的句柄写法

1. f= @(x) 2*(1600 -40*r*cos(x))./(r^2-80*r*cos(x)+1600+z^2).^1.5;

复制代码

应该是

1. 2*(1600-40*40*cos(x))/(40^2-80*40*cos(x)+1600+（-4.5）^2)^1.5

复制代码

所以只需要把原来的式子改为

1. f=strcat('2*(1600-40*',num2str(r),'*cos(x))/(',num2str(r),'^2-80*',num2str(r),'*cos(x)+1600+(',num2str(z),')^2)^1.5');
   
2. f = inline(vectorize(f));

复制代码

就没有问题了

qibbxxt

另外我觉得你的程序，可以写的简洁一些，当然是我的个人看法

1. rr = linspace(40.00 , 45.00 , 100) ;          %产生r=40到45的数组
   
2. zz = linspace(-4.5 , 4.5 , 6) ;               %产生每匝的z位移从-4.5到4.5
   
3. f = @(r,z)quad(@(x) 2*(1600 -40*r*cos(x))./(r^2-80*r*cos(x)+1600+z^2).^1.5,0,pi);
   
4. [rz{[2 1]}]=meshgrid(zz,rr);
   
5. b = arrayfun(f,rz{:});
   
6. b=[b,sum(b,2),rr'];

复制代码