不等距离散点的积分程序

sigma665

1. function [result error ind]= cubint(X,F,N,ia,ib)
   
2. % 参考资料：Philio J. Davis 数值积分法 高等教育出版社 1986 p296-298
   
3. % X,F:积分点，函数值
   
4. % N: 积分点数，N>=4
   
5. % ia ib: ia,ib须为（1，N）整数 X(ia),X(ib)
   
6. % result: 积分值
   
7. % error: 误差估计
   
8. % ind: 条件满足为1，否则为0
   
9. % 对离散点三次插值，给出了基于四阶差商逼近四阶导数的误差估计
   
10. %
    
11. % 例：
    
12. % x=[0,0.1,0.15,0.2,0.23,0.3,0.4,0.45,0.48,0.53,0.62,0.78,0.82,0.89,0.92,1];
    
13. % y=x.^0.25;
    
14. % [result error ind]= cubint(x,y,16,1,16)；
    
15. % 运行结果：
    
16. % result = 0.7905 精确解：0.8
    
17. % error = -8.7208e-005
    
18. % ind = 1
    
19. 
    
20. 
    
21. ind=0;
    
22. 
    
23. if N<4 || ia<1 || ib>N
    
24. disp('wrong in N ia ib');
    
25. else
    
26. ind=1;
    
27. result=0;
    
28. error=0;
    
29. 
    
30. if ia==ib
    
31. disp('ia=ib,please check ia or ib');
    
32. else
    
33. if ia>ib
    
34. ind=-1;
    
35. it=ib;
    
36. ib=ia;
    
37. ia=it;
    
38. end
    
39. s=0;
    
40. c=0;
    
41. r4=0;
    
42. J=N-2;
    
43. if ia<(N-1) || N==4
    
44. J=max(3,ia);
    
45. end
    
46. K=4;
    
47. if ib>2 || N==4
    
48. K=min(N,ib+2)-1;
    
49. end
    
50. for i=J:K
    
51. if i<=J
    
52. H2=X(J-1)-X(J-2);
    
53. d3=(F(J-1)-F(J-2))/H2;
    
54. H3=X(J)-X(J-1);
    
55. d1=(F(J)-F(J-1))/H3;
    
56. H1=H2+H3;
    
57. d2=(d1-d3)/H1;
    
58. H4=X(J+1)-X(J);
    
59. r1=(F(J+1)-F(J))/H4;
    
60. r2=(r1-d1)/(H4+H3);
    
61. H1=H1+H4;
    
62. r3=(r2-d2)/H1;
    
63. if ia<=1
    
64. result=H2*(F(1)+H2*(0.5*d3-H2*(d2/6-(H2+H3+H3)*r3/12)));
    
65. s=-H2^3*(H2*(3*H2+5*H4)+10*H3*H1)/60;
    
66. end
    
67. else
    
68. H4=X(i+1)-X(i);
    
69. r1=(F(i+1)-F(i))/H4;
    
70. r4=H4+H3;
    
71. r2=(r1-d1)/r4;
    
72. r4=r4+H2;
    
73. r3=(r2-d2)/r4;
    
74. r4=(r3-d3)/(r4+H1);
    
75. end
    
76. if i>ib || i<=ia
    
77. error=error+r4*s;
    
78. else
    
79. result=result+H3*((F(i)+F(i-1))*0.5-H3*H3*(d2+r2+(H2-H4)*r3)/12);
    
80. c=H3^3*(2*H3^2+5*(H3*(H4+H2)+2*H2*H4))/120;
    
81. error=error+(c+s)*r4;
    
82. if i==j
    
83. s=s+c+c;
    
84. else
    
85. s=c;
    
86. end
    
87. end
    
88. if i>=K
    
89. if ib>=N
    
90. result=result+H4*(F(N)-H4*(0.5*r1+H4*(r2/6+(H3+H3+H4)*r3/12)));
    
91. error=error-H4^3*r4*(H4*(3*H4+5*H2)+10*H3*(H2+H3+H4))/60;
    
92. else
    
93. if ib>=N-1
    
94. error=error+s*r4;
    
95. end
    
96. end
    
97. else
    
98. H1=H2;
    
99. H2=H3;
    
100. H3=H4;
     
101. d1=r1;
     
102. d2=r2;
     
103. d3=r3;
     
104. end
     
105. end
     
106. if ind~=1
     
107. it=ib;
     
108. ib=ia;
     
109. ia=it;
     
110. result=-result;
     
111. error=-error;
     
112. ind=1;
     
113. end
     
114. end
     
115. end
     
116. end

复制代码

最近在做边界元,要积分
也看到有人问过离散点的积分
刚好,刚调完,结果正确
是根据fortran改过来的
抛砖引玉了

[ 本帖最后由 sigma665 于 2008-1-16 20:19 编辑 ]