数值微分和数值积分
数值微分
(1)数值差分与差商
可以用差商近似计算导数值:
(2)数值微分的实现
diff() 向前差分函数
dx=diff(x):
计算向量x的一阶向前差分,dx[i]=x[i+1]-x[i],i=1,2,3…n-1.dx=diff(x,n):
计算向量x的n阶向前差分。dx=diff(A,n,dim):
计算矩阵A的n阶差分。dim=1时按列计算差分,=2时按行计算差分。
举个例子
x=[0,sort(2*pi*rand(1,5000)),2*pi];
y=sin(x);
f1=diff(y)./diff(x)%利用差商来近似导数
f2=cos(x(1:end-1))%实际导数
plot(x(1:end-1),f1,x(1:end-1),f2,'r:')
两者较接近
数值积分
原理
函数调用
[l,n]=quad(filename,a,b,tol,trace)
基于自适应辛普森方法[l,n]=quadl(filename,a,b,tol,trace)
基于自适应GL方法(效率比较高)
参数含义:I=integral(filename,a,b)
基于全局自适应积分方法。I是计算得到的积分。[I,err]=quadgk(filename,a,b)
基于自适应高斯—克朗罗得方法。err返回近似误差。I=trapz(x,y)
x自变量,y因变量。
举个例子
format long
f=@(x) 4./(1+x.^2);
[I,n]=quad(f,0,1,1e-8)
[I,n]=quadl(f,0,1,1e-8)
(atan(1)-atan(0))*4
结果
I =3.141592653733437
n =61
I =3.141592653589806
n =48
ans =3.141592653589793
f=@(x) 1./(x.*sqrt(1-log(x).^2));
I=integral(f,1,exp(1))
结果:
I =1.570796326795570
f=@(x) sin(1./x)./x.^2;
I=quadgk(f,2/pi,+Inf)
结果:
I = 1.000000000000000
x=1:6;
y=[6 8 11 7 5 2];
plot(x,y,'-ko');
grid on
axis([1 6 0 11]);
I1=trapz(x,y)
I2=sum(diff(x).*(y(1:end-1)+y(2:end))/2)
结果:
I1 =35
I2 =35
下图是x—y图像