一种积分限以及被积函数中含有未知数的积分方程的MATLAB解法

rocwoods

(发在研学上的帖子)
最近帮助朋友解决一个问题，遇到了求解积分方程。在很多论坛上翻了翻帖子发现对这个问题的解决大多停留在符号求解或者转化为微分方程再求解上。由于很多积分没有解析解，因此符号求解不仅费时而且往往要碰壁。（我曾经符号求解过80多个积分方程，可能是比较复杂，结果程序一晚上才解了8个左右，速度令人无法忍受）
而转化为微分方程需要根据方程的形式（有的还不好化），且需要一定的数学基础，通用性较差。
为此，给出一种可行的数值解法，希望抛砖引玉，大家一起完善总结积分方程的解决办法。
这种方法主要利用MATLAB中的fzero函数以及匿名函数。fzero函数所用的算法是zeroin算法，该算法非常简单而且安全，其核心思想是将二分法的可靠性与逆二次插值法、割线法结合起来。这个算法是比较简单而且安全的，非常适合求在零点变号的函数的零点。另外利用匿名函数（7.0以上版本支持）可以方便的建立被积函数以及积分方程。
譬如针对如图所示的积分方程：
代码如下：
tt=@(s) ['(x-' num2str(s) ').*exp(-x.^2/2)'];
ff=@(s) quadl(tt(s),s,10)-1/2;
sol=fzero(ff,3)
求得的结果为sol =
0.49458500000000

frensel

匿名函数太慢，最好写成nested function的形式，即：

1. 
   
2. % rocwoods.m
   
3. function y = rocwoods(s)
   
4. y = quadl(@fun, s, 10) - 1 / 2;
   
5.     function f = fun(x)
   
6.         f = (x - s) .* exp(-x.^2 / 2);
   
7.     end
   
8. end
   

复制代码

1. 
   
2. >> fzero(@rocwoods, 3)
   
3. ans =
   
4.    0.494589233651044
   

复制代码

rocwoods

呵呵，不错。我只知道匿名函数效率比内联函数高，对嵌套函数和匿名函数的效率没有比较过。
去你的空间看了看，发现你对分布式计算挺有研究，以后遇到这方面的问题还要多指教^_^

frensel

原帖由 rocwoods 于 2007-3-23 09:47 发表
呵呵，不错。我只知道匿名函数效率比内联函数高，对嵌套函数和匿名函数的效率没有比较过。
去你的空间看了看，发现你对分布式计算挺有研究，以后遇到这方面的问题还要多指教^_^

其实nested function也很慢。我有些工作要调用同一个积分上万次，这个积分又包含其他的未知数，使函数不得不写成nested function的形式。最后测试，程序大部份时间都耗在这函数的调用上，改写成匿名函数效率又更低，差距非常明显。

rocwoods

不知道有无更高效率的办法。期待潜水的各位大牛完善这个问题。

frensel

原帖由 rocwoods 于 2007-3-23 14:25 发表
不知道有无更高效率的办法。期待潜水的各位大牛完善这个问题。

这个问题我在这发过贴问，但没解决方法……
可看这：http://www.simwe.com/forum/viewthread.php?tid=764661&highlight=%2Bfrensel

dty7611477

(发在研学上的帖子)
最近帮助朋友解决一个问题，遇到了求解积分方程。在很多论坛上翻了翻帖子发现对这个问题的解决大多停留在符号求解或者转化为微分方程再求解上。由于很多积分没有解析解，因此符号求解不仅费时而 ...
rocwoods 发表于 2007-3-22 10:52

积分限以及被积函数中含有未知数的积分,対上式求最小值优化，优化变量为a b c
望高人指点！

dty7611477

7# dty7611477
function y=opt_fun(x)
%数值积分
y=quadl(@t,xia(x),shang(x));
%求积分上限
function th_x=shang(x)
s=x(3)+x(4);
th3=-asin((x(3).^2-x(4).^2+s.^2)./(2*s.*x(3)));
g=-2*x(1).*(x(6)+x(3).*sin(th3));
h=-2*x(1).*(x(5)+x(3).*cos(th3));
k=x(3).^2+x(1).^2+2*x(3).*x(5).*cos(th3)+2*x(3).*x(6).*sin(th3)+x(5).^2+x(6).^2-x(1).^2;
th_x=2*atan((sqrt(d^2+e^2-f^2)+d)/(e-f));
end
%求积分下限
function th_6=xia(x)
s=x(3)+x(4)-0.006;
th3=-asin((x(3).^2-x(4).^2+s.^2)./(2*s.*x(3)));
gg=-2*x(1).*(x(6)+x(3).*sin(th3));
hh=-2*x(1).*(x(5)+x(3).*cos(th3));
kk=x(3).^2+x(1).^2+2*x(3).*x(5).*cos(th3)+2*x(3).*x(6).*sin(th3)+x(5).^2+x(6).^2-x(1).^2;
th_6=2*atan((sqrt(d^2+e^2-f^2)+d)/(e-f));
end
%被积函数
function T=t(th)
d=2*x(2).*(x(1).*sin(th)-x(6));
e=2*x(2).*(x(1).*cos(th)+x(5));
f=(x(1).*cos(th)+x(5)).^2+(x(1).*sin(th)-x(6)).^2+x(2).^2-x(3).^2;
th2=2*atan((sqrt(d.^2+e.^2-f.^2)+d)./(e-f));

a=2*x(3).*(x(6)-x(1).*sin(th));
b=-2*x(3).*(x(5)+x(1).*cos(th));
c=x(3).^2-x(2).^2+(-x(1).*sin(th)).^2+(x(5)+x(1).*cos(th)).^2;
th3=2*atan((sqrt(a.^2+b.^2-c.^2)+d)./(b-c));
th4=acos(-x(3).*cos(th3)./x(4));
%T是最终的被积函数，上面是th2，th3，th4关于th的函数，要对th积分
T=x(1).*cos(th)+(x(1).*sin(th3-th).*sin(th4).*cos(th)-x(1).*sin(th2-th).*sin(th3).*sin(th4))./(sin(th2-th3).*sin(th4));
end

end
/////////////////////////////////////////////////////////////////////////////////////////////////////////
%利用工具箱中的函数fmincon优化，求最优的L1 L2 L3 L4 x y使T最小值
x0=[0.06;0.347;0.114;0.210;-0.308; 0.090];
[x,fim] = fmincon(@opt_fun,x0,[],[],[],[],[],[],@constrion)
  /////////////////////////////////////////////////////////////////////////////////////////////////////////
%约束函数
function [c1,c2]=constrain(x)
c1=[x(1)+sqrt(x(5)^2+x(6)^2)-x(2)-x(3);x(1)-sqrt(x(5)^2+x(6)^2)+x(2)-x(3);x(1)-sqrt(x(5)^2+x(6)^2)-x(2)+x(3)];
c2=[];
这个积分和1楼的很相似，就是上下限都是很长的式子，我用函数表示的，不知道哪里出错，望高手指点，谢谢

waynebuaa

对于这类似的问题，最好先自己化简一下，换成易于处理的形式。 比如，楼主的积分，可以借用误差函数erf(x)，把积分过程去掉，再来求根
                        

TBE_Legend

A general insight of computational method of this kind of problem can be found in many papers dealing with solving  vloterra integro-differential equations. it seems not too difficult.

waynebuaa

volterra integro-differential equations?

sounds an old jargon , thank you