资源描述
计算向量、矩阵差分
> x=1:5
> x
> diff(x)
[1] 1 1 1 1
> y=1:12
> z=matrix(y,3,4,T)
> z
> diff(z)
[,1] [,2] [,3] [,4]
[1,] 4 4 4 4
[2,] 4 4 4 4
> diff(x,diff=2)
[1] 0 0 0
#求一元函数导数
f=expression((5*x)/(1+x^2))
D(f,"x")
f=expression(2*x^2)
D(f,"x")
#求二元函数偏导数及梯度
f=deriv(z~x^2+y^2+x*y,c('x','y'),func=T, hessian=F)
f
f(1,2)
例1梯形积分
trape=function(fun,a,b,tol=1e-6){
N=1;h=b-a
T=h/2*(fun(a)+fun(b))
repeat{
h=h/2;x=a+(2*1:N-1)*h
I=T/2+h*sum(fun(x))
if(abs(I-T)<tol) break
N=2*N;T=I}
I}
source("trape.R")
f=function(x) exp(-x^2)
trape(f,-1,1)
高精度积分
f=function(x) exp(-x^2)
integrate(f,-1,1)
例2简单二重积分
library(cubature)
f<-function(x){
1-x[1]/4-x[2]/3
}
adaptIntegrate(f, c(-2, -1), c(2, 1), maxEval=10000)
例3
直接求和计算
n=5;h=seq(0.1,0.9,by=0.2)
fun=function(h)(2-(exp(2*h)+exp(-2*h))/5)
r=fun(h)
r
w=pi*1/n*sum(r^2)
w
转化为定积分
f=function(h)(2-(exp(2*h)+exp(-2*h))/5)^2
I=pi*integrate(f,0,1)
as.character(I)
I=as.numeric(I[1])
w=I*pi
w
例4
R<-5
b<-0.1
g<-9.81
z1<-0.5-R
z2<-R
n<-100
h<-(z2-z1)/n
#z=z1:h:z2
f<-function(z)(R^2-z^2)/(b^2*sqrt(2*g*(z+R)))
I=integrate(f,z1,z2)
I
as.character(I)
I=as.numeric(I[1])
t=I/60/60
t
展开阅读全文