这个积分方程matlab应该怎么求？
xco与xc1是两个未知数的函数
方程左边的积分部分尝试用quadl函数来表达
fun1=inline('k.*sqrt((1.521145+x(1).*(1-erf(x/(2.*(sqrt(x(2).*t))))))^2-N1^2)');
R1=quadl(fun1,0,x0,1e-8);
fun2=inline('k.*sqrt((1.521145+x(1).*(1-erf(x/(2.*(sqrt(x(2).*t))))))^2-N2^2)');
R2=quadl(fun2,0,x1,1e-8);
但是运行出错
??? Error using ==> inline.feval
Not enough inputs to inline function.
Error in ==> quadl at 64
y = feval(f,x,varargin{:}); y = y(:).';
Error in ==> only at 14
R1=quadl(fun1,0,x0,1e-8);
Error in ==> fsolve at 180
fuser = feval(funfcn{3},x,varargin{:});
哪位能帮忙看看是怎么回事吗？
去掉inline函数 改用下面的程序
R1=quadl('k*sqrt((1.521145+x(1)*(1-erf(x/(2*(sqrt(x(2)*t))))))^2-N1^2)',0,x0,1e-8);
R2=quadl('k*sqrt((1.521145+x(1)*(1-erf(x/(2*(sqrt(x(2)*t))))))^2-N2^2)',0,x1,1e-8);
运行结果一样
[ 本帖最后由 ChaChing 于 2009-6-21 15:48 编辑 ] x0给什麽? x0在主函数fsolve函数里边给了初值的
'k*sqrt((1.521145+x(1)*(1-erf(x/(2*(sqrt(x(2)*t))))))^2-N1^2)'这个式子里边含有erf（）函数
程序调试到这一步就停了，是不是对erf函数没法直接积分呢？那要怎么改
******only.m函数******
function y=only(x)
t=1200; N0=1.554336; N1=1.536857; ns=1.521145;
k=2*pi/(632.8e-9); B0=k*N0; B1=k*N1;
x0=2*finverse_erfc((N0-ns)/x(1))*sqrt(x(2)*t);
x1=2*finverse_erfc((N1-ns)/x(1))*sqrt(x(2)*t);
ierfcx1=int('k*sqrt((1.521145+x(1)*erfc(z/(2*(sqrt(x(2)*t)))))^2-N0^2)','z',0,x0);
ierfcx2=int('k*sqrt((1.521145+x(1)*erfc(z/(2*(sqrt(x(2)*t)))))^2-N1^2)','z',0,x1);
R1=subs(ierfcx1); R2=subs(ierfcx2);
tg1=sqrt((B0^2-k^2)/(k^2*(ns+x(1))^2-B0^2))-x(1)*(k^2)*(ns+x(1))/(2*sqrt(pi*x(2)*t)*(k^2*(ns+x(1))^2-B0^2)^1.5);
tg2=sqrt((B1^2-k^2)/(k^2*(ns+x(1))^2-B1^2))-x(1)*(k^2)*(ns+x(1))/(2*sqrt(pi*x(2)*t)*(k^2*(ns+x(1))^2-B1^2)^1.5);
eq1=R1-pi/4-atan(tg1); eq2=R2-5*pi/4-atan(tg2);
y=；
**********finverse_erfc.m函数********
function c=finverse_erfc(z);
% count inverse function of erfc
% author email of the program:% zjliu2001@163.com
a=600; b=0;
if z>1 | z<0;
error([[' z takes a error value!'],...
['; please check it''s value!']]);
end
c=(a+b)/2;
while abs(erfc(c)-z)>1e-6;
c=(a+b)/2;
if erfc(c)-z<0, a=c; else b=c; end
end
********主函数******
x0=; options=optimset('display','iter');
=fsolve('only',x0,options)
[ 本帖最后由 ChaChing 于 2009-6-21 15:54 编辑 ] 楼主还用了萝卜写的程序了啊？呵呵，你写的程序得需要较大的改动，晚上有时间可以帮你想想这个问题。 楼上说的是finverse_erfc.m函数啊？
呵呵，在网上看到的就下下来了，感觉还挺好的！
嗯，谢谢recwoods~
调试时发现，运行到这一步的时候
ierfcx1=int('k*sqrt((1.521145+x(1)*erfc(z/(2*(sqrt(x(2)*t)))))^2-N0^2)','z',0,x0);
ierfcx2=int('k*sqrt((1.521145+x(1)*erfc(z/(2*(sqrt(x(2)*t)))))^2-N1^2)','z',0,x1);
fsolve设定的x初值没有传递进int积分函数
这个是为什么？应该怎么改呢？
[ 本帖最后由 ChaChing 于 2009-6-21 15:56 编辑 ] 建议楼主only函数用nested function形式写，这样可以解决参数传递问题。楼主可以查看帮助文档或者google、baidu了解nested function用法。
关于向被积函数传参数，楼主可以参考这个帖子。
http://forum.vibunion.com/forum/thread-42369-1-1.html 谢谢rocwoods~~
研究了下您说的nested function，
有个问题想问一下
我这里边需要编两个嵌套函数，s1跟s2是同级的。
那我可以这样使用吗？
function r = fatherfunction()
blah
function s1 = childfunction()
blah
function s2 = childfunction()
end
end
end function f1=f1(z)
f1=sqrt(k^2.*(ns+x(1).*erfc(z/(2*(x(2)*t)^0.5))).^2-B0^2);
end
function f2=f2(z)
f2=sqrt(k^2.*(ns+x(1).*erfc(z/(2*(x(2)*t)^0.5))).^2-B1^2);
end
ierfcx1=quadl(@f1,0,x0); ierfcx2=quadl(@f2,0,x1);
程序运行到f1=sqrt(k^2.*(ns+x(1).*erfc(z/(2*(x(2)*t)^0.5))).^2-B0^2);就停住了
x0=; options=optimset('display','iter');
=fsolve('only',x0,options)
Norm of First-order Trust-region
IterationFunc-count f(x) step optimality radius
0 3 2789.92 2.49e+013 1
MATLAB一直显示busy。。。
[ 本帖最后由 ChaChing 于 2009-6-21 16:00 编辑 ]
回复 7楼 悠若谷 的帖子
你的blah是什么?在另一个帖子里如何求erfc(x)函数的平方的积分我已经给你回了，你可以参考下。你可以把你卡住的问题以最简单的例子呈现出来，不是光学专业的，不懂你那个方程各个物理量含义，方程看起来就费劲，不过我觉得你卡住的地方是可以用简单的例子说明白的。还有，工作后很少用QQ了，有问题可以在论坛上提，这样讨论的结果可以方便记录下来，同时方便以后可能遇到类似问题的人，这也是论坛的初衷之一。[ 本帖最后由 ChaChing 于 2009-6-21 16:03 编辑 ] 改了一下现在不是卡住了，设置断点单步运行的时候fsolve设定的初值传进去没有问题，第一,二次搜索很正常，但是到第三次搜索的时候程序就中断了。。。我了解的matlab知识不够用的，不知道是不是方程太复杂了不适合用fsolve。
Norm of First-order Trust-region
IterationFunc-count f(x) step optimality radius
0 3 9.66525 2.03e+012 1
??? Input must be real.
Error in ==> erfc at 18
y = erfcore(x,1);
Error in ==> only>f1 at 14
f1=sqrt(k^2.*(ns+x(1)*erfc(z/(2*(x(2)*t)^0.5))).^2-B0^2);
Error in ==> quad at 62
y = f(x, varargin{:});
Error in ==> only at 16
R1=quad(@f1,0,x0);
Error in ==> optim\private\trustnleqn at 210
F = feval(funfcn{3},x,varargin{:});
Error in ==> fsolve at 295
=...
Error in ==> fmain at 3
=fsolve('only',x0,options)
[ 本帖最后由 悠若谷 于 2009-6-12 14:29 编辑 ] 这个程序，如果把复杂的地方拿出来用简单的例子表示基本上运行是没问题的，但是都凑到一起就不行了。您要是有空的话帮我看看吧，这个程序的思想很简单，就是变量很多。
%******fmain函数*****
x0=;
options=optimset('display','iter');
=fsolve('only',x0,options)
%*****only.m********
function y=only(x)
t=1200;
N0=1.554336;
N1=1.536857;
ns=1.521145;
k=2*pi/(632.8e-9);
B0=k*N0;
B1=k*N1;
x0=2*finverse_erfc((N0-ns)/x(1))*sqrt(x(2)*t);
x1=2*finverse_erfc((N1-ns)/x(1))*sqrt(x(2)*t);
tg1=sqrt((B0^2-k^2)/(k^2*(ns+x(1))^2-B0^2))-x(1)*(k^2)*(ns+x(1))/(2*sqrt(pi*x(2)*t)*(k^2*(ns+x(1))^2-B0^2)^1.5);
tg2=sqrt((B1^2-k^2)/(k^2*(ns+x(1))^2-B1^2))-x(1)*(k^2)*(ns+x(1))/(2*sqrt(pi*x(2)*t)*(k^2*(ns+x(1))^2-B1^2)^1.5);
function f1=f1(z)
f1=sqrt(k^2.*(ns+x(1)*erfc(z/(2*(x(2)*t)^0.5))).^2-B0^2);
end
R1=quad(@f1,0,x0);
function f2=f2(z)
f2=sqrt(k^2.*(ns+x(1)*erfc(z/(2*(x(2)*t)^0.5))).^2-B1^2);
end
R2=quad(@f2,0,x1);
eq1=R1-pi/4-atan(tg1);
eq2=R2-5*pi/4-atan(tg2);
y=;
end
%***finverse_erfc.m******
function c=finverse_erfc(z);
% count inverse function of erfc
% author email of the program:% zjliu2001@163.com
a=600;
b=0;
if z>1 | z<0;
error([[' z takes a error value!'],...
['; please check it''s value!']]);
end
c=(a+b)/2;
while abs(erfc(c)-z)>1e-6;
c=(a+b)/2;
if erfc(c)-z<0;
a=c;
else
b=c;
end
end
页:
[1]