带2参数的数值积分
本帖最后由 yxlnbu 于 2011-3-10 10:18 编辑大家好，我有一个函数'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
functionb = 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为你要的数据 ！ 回复 2 # zhouyang664 的帖子
谢谢您，在这个问题上花了好长时间不得其解，谢谢您的帮助。能解释下为什么inline函数不行么？有时候就说输入不够之类的。 回复 3 # yxlnbu 的帖子
LZ比較下兩者差異, 應該即可體會
f='2*x+3';
ff=inline('f')
ff=inline(f) 回复 4 # ChaChing 的帖子
谢谢，刚才我试了下，一个是以x为变量，另一个是以f为变量的，请问是这样理解么？另外，我把上面的程序去掉引号后，程序计算就出现了奇异，而应用zhouyang664的程序的时候就没有出现这个问题。但把rr和zz值缩小十倍之后也出现了奇异
Warning: Maximum function count exceeded; singularity likely. yxlnbu 发表于 2011-3-11 10:32 static/image/common/back.gif
回复 4 # ChaChing 的帖子
谢谢，刚才我试了下，一个是以x为变量，另一个是以f为变量的，请问是这样理解么 ...
没看懂你的问题
直接把你修改后而出现问题的代码贴上来，是查找问题的最直接方法 回复 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
楼主可以自己用断点调试查一下问题的所在啊！ 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. 回复 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) ;
软件的警讯是不会错的! 也就是说函数有奇异, 但会因使用inline或@而不同吗?
我看楼主也是在休周末了, 就留下问题, 给有兴趣自我提升的人试试吧! 回复 11 # ChaChing 的帖子
谢谢各位，我是回家了几天，没有上论坛了。今天刚回来。
函数中有改变的就是@和inline函数调用，还有的就是‘num2str’ 函数的调用。大家可以试试。 本帖最后由 ChaChing 于 2011-3-15 21:21 编辑
都忘了这了! 我想没人想试看看赚威望了
看看差异!!
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');
回复 13 # ChaChing 的帖子
按说Matlab不会出现这种错误的，其实这两种写法，函数是不一样的 f=strcat('2*(1600-40*',num2str(r),'*cos(x))/(',num2str(r),'^2-80*',num2str(r),'*cos(x)+1600+',num2str(z),'^2)^1.5');这种写法，当z是负值的时候，就是 2*(1600-40*40*cos(x))/(40^2-80*40*cos(x)+1600+-4.5^2)^1.5而后面的句柄写法 f= @(x) 2*(1600 -40*r*cos(x))./(r^2-80*r*cos(x)+1600+z^2).^1.5;应该是 2*(1600-40*40*cos(x))/(40^2-80*40*cos(x)+1600+（-4.5）^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');
f = inline(vectorize(f));就没有问题了 另外我觉得你的程序，可以写的简洁一些，当然是我的个人看法rr = linspace(40.00 , 45.00 , 100) ; %产生r=40到45的数组
zz = linspace(-4.5 , 4.5 , 6) ; %产生每匝的z位移从-4.5到4.5
f = @(r,z)quad(@(x) 2*(1600 -40*r*cos(x))./(r^2-80*r*cos(x)+1600+z^2).^1.5,0,pi);
}]=meshgrid(zz,rr);
b = arrayfun(f,rz{:});
b=;
页:
[1]
2