Simpson 法求数值微分

１、编程实现以下科学计算算法，并举一例应用之。（参考书籍《精 通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９ 年）

“Simpson 法求数值微分” 算法说明

辛普森数值微分是用来求等距节点在节点处的导数的，辛普森数值微分公式如下： 4 1 f’(x1) 3(y2-y0)/h-f’(x0) 1 4 1 f’(x2) 3(y3-y1)/h 1 4 . 3(y4-y2)/h f’(x3)

. . . . . . . . . . . . . . . . 1 . .

1 4 3(yn-yn-2)/h-f’(xn)

f’(x=n-1)

其中，yn=f(xn),xn=xn+nh

如果端点导数值－f(x0)和-f’(xn)未知，则将它们的中点微分公式近似，这时的辛普森数值微分公式为：

2 0 f’(x1) y2-y0/h

1 4 1 f’(x2) 3(y3-y1)/h 1 4 . 3(y4-y2)/h f’(x3) . . . . . . . . . . . .

. . . . 1 . . 0 .2 f’(x=n-1) (yn-yn-2)/h

在MATLAB中，编程实现的辛普森数值法（用于已知函数表达式）函数为：CISimpson。 功能：辛普森数值法求已知函数在某点的导数值。 调用格式：df=CISimpson(func,x0,n,h)。 其中，func为函数名； xo为求导点；

n为将已知函数离散的数据点数；

h为离散步长； df为导数值。

辛普森数值微分法（适用于已知函数表达式）的MATLAB程序代码如下： function df=CISimpson(func,x0,n,h)

%辛普森数值法求已知函数func在x0点的导数值 %函数名：func %求导点：x0

%将已知函数离散的数据点数：n %离散步长：h %导数值：df

= = if nargin ==2 %以下是参数的判断过程 h=0.1; n=5; else

if(nargin ==3) if (n<5)

disp('n不能小于5!'); return; else

h=0.1; end

else (nargin ==4 && h ==0.0) disp('h不能为0!'); return; end end

for(i=1:n) %这是个循环计算节点的函数值 if (mod(n,2) ==0)

y(i)=subs(sym(func),findsym(sym(func)),x0+(i-n/2)*h); else

y(i)=subs(sym(func),findsym(sym(func)),x0+(i-(n+1)/2)*h); end end

f(1)=(y(3)-y(1))/(2*h);

f(2)=(y(n)-y(n-2))/(2*h); %这两行用中心微分法给出端点的导数 b(1:n-2,1)=zeros(n-2,1); b(1,1)=3*(y(3)-y(1))/h-f(1); b(n-2,1)=3*(y(n)-y(n-2))/h-f(2); for (i=2:(n-3))

b(i,1)=3*(y(i+2)-y(i))/h;

end %这一块是辛普森公式的右边的列向量 for(i=1:n-2) for(j=1:n-2)

if((i ==j+1) ||(j ==i+1)) A(i,j)=1; else if(i ==j) A(i,j)=4; end end end

end %这一块是系数矩阵 [Q,R]=qr(A);

DF = R\\(Q\\b); %用QR分解法求解 if(mod(n,2) ==0) df = DF(n/2); else

df=DF((n+1)/2);

end %这里是求出x0处的导数值

辛普森数值微分法（适用于数据向量）的另一个MATLAB程序代码如下： function df=DISimpson(X,Y,n,p)

%辛普森数值法求n个数据点在第p个点处的导数值 %离散数据的x坐标向量：X %离散数据的y坐标向量：Y %数据点数：n %要求导数的点:p %导数值:df if n<5

disp('n不能小于5!'); return; end

if p==0

disp('n不能小于0!'); return; end

h=X(2)-X(1);

xx=linspace(X(1),X(n),h); if(xx~=X)

disp('节点之间不是等距的！'); return; end

f(1)=(Y(3)-Y(1))/(2*h); f(2)=(Y(n)-Y(n-2))/(2*h); b(1,1)=3*(Y(3)-Y(1))/h-f(1); b(n-2,1)=3*(Y(n)-Y(n-2))/h-f(2); for(i=2:(n-3))

b(i,1)=3*(Y(i+2)-Y(i))/h; end

for(i=1:n-2) for(j=1:n-2)

if((i==j+1)||(j==i+1)) A(i,j)=1; else if(i==j) A(i,j)=4; end end end end

[Q,R]=qr(A); DF=R\\(Q\\b); if(p==1) df=f(1); else

df=DF(p-1);

end

辛普森数值微分法应用举例

分别采用5,10,100个离散点的辛普森法求函数exp（x）在x=2.5的导数值 运算结果