学习R语言编程——常用算法——导数与微积分的近似计算

#########################常用算法——导数与微分的近似计算#######################

###例题：用导数的近似计算求函数f(x) = -4*x^2+3*x+2在[0,1]的极值

f = function(x) -4*x^2+3*x+2

ub = 1

lb = 0

m = 1000000

step = (ub-lb)/m

x = 0

for(i in 1:1000000){

df = (f(x+step)-f(x))/step

if(abs(df)<1e-5){cat("x is ",x,"extreme value is ",f(x))}

x = x+step

}

#矩阵方法

f = function(x) -4*x^2+3*x+2

x = seq(0,1,by = .00001)

x0 = x[which.max(f(x))]

x0;f(x0)

 

###用导数的近似计算求出函数f(x) = (1/3)*x^3-(1/2)*x^2-2*x在[-2,3]的所有极值点。

#判定方法1

f = function(x) (1/3)*x^3-(1/2)*x^2-2*x

x = -2

step = 1/10000

for(i in 1:50000){

fd = (f(x-step)-f(x))*(f(x+step)-f(x))  #判断是否是极值点

if(fd>0)print(x)

x = x+step

}

#判定方法2——定义判定

f = function(x) (1/3)*x^3-(1/2)*x^2-2*x

ub = 3

lb = -2

m = 5000000

step = (ub-lb)/m

x = -2

for(i in 1:m){

df = (f(x+step)-f(x))/step

if(abs(df)<1e-5)cat("extreme value point is ",x,'\n')

x = x+step

}

 

 

 

########################常用算法——定积分的近似计算###########################

###利用矩形法近似计算积分shit(0,3)x^2dx

#利用循环结构

f = function(x) x^2

x = 0

lb = 0

ub = 3

m = 30000

step = (ub-lb)/m

s = 0

for(i in 1:30000){

s = s + step*f(x)

x = x + step

}

cat(s)

#矩阵方法

f = function(x) x^2

x = seq(0,3,length = 30001)

a = x[2:length(x)]

sum(f(a)*(3/length(a)))                     #矩形法

 

###求积分shit(0,4)(x+2)/sqrt(2*x+1)dx

#矩阵方法

f = function(x) (x+2)/sqrt(2*x+1)

x = seq(0,4,by = .001)

a = x[2:(length(x)-1)]

sum(f(a)*.001)                              #矩形法

#方法2

a=seq(0,4,by=1e-3)

tail(a)

length(a)

f=function(x) (x+2)/sqrt(2*x+1)

sum(f(a[-length(a)])*1e-3)                  #矩形法

 

###利用梯形法近似计算积分shit(0,1)exp(-1*x^2)dx

#利用循环结构求解

f = function(x) exp(-1*x^2)

lb = 0

ub = 1

m = 1000

step = (ub - lb)/m

x = 0

s = 0

for(i in 1:m){

long = (f(x)+f(x+step))/2

s = s + step*long

x = x+step###注意这里x也要更新，一步一步递增前进

}

cat(s)

 

#矩阵方法

f = function(x) exp(-1*x^2)

x = seq(0,1, by = .001)

a = x[1:(length(x)-1)]

b = x[2:length(x)]

c = (a+b)/2

s = sum(f(c)*.001)###注意这里是f(c)即函数值与步长相乘

s

 

###用矩形法和梯形法近似计算积分shit(0,1)1/(1+x^2)dx,

###并由此计算pi的近似值，n = 6被积函数值取5为小数

#矩形法——循环结构求解

f = function(x) 1/(1+x^2)

lb = 0

ub = 1

m = 1000

step = (ub - lb)/m

x = 0

s = 0

for(i in 1:m){

s = s+f(x)*step

x = x+step

}

p = 4*s;p

round(p,digits = 5)

 

#矩形法——矩阵方法求解

f = function(x) 1/(1+x^2)

x = seq(0,1,by = .001)

a = x[1:(length(x)-1)]

s = sum(f(a)*step);s

p = 4*s;p

round(p,digits = 5)

 

#梯形法——循环结构求解

f = function(x) 1/(1+x^2)

lb = 0

ub = 1

m = 1000

step = (ub - lb)/m

x = 0

s = 0

for(i in 1:m){

s = s + ((f(x)+f(x+step))/2)*step

x = x+step

}

p = 4*s;p

round(p,digits = 5)

 

#梯形法——矩阵方法

f = function(x) 1/(1+x^2)

x = seq(0,1,by = .001)

a = x[1:(length(x)-1)]

b = x[2:length(x)]

c = (a+b)/2

s = sum(f(c)*step)

p = 4*s

round(p,digits = 5)

 

###定积分在几何上的应用——求弧长

#使用循环结构求解

f = function(x) log((1+sqrt(1-x^2))/x)-sqrt(1-x^2)

lb = 1/2

ub = 1

m = 10000

arc = 0

step = (ub-lb)/m

x = 1/2

for(i in 1:m){

arc = arc + sqrt(step^2 + (f(x+step)-f(x))^2)

x = x + step

}

cat(arc)

#使用矩阵方法求解

f = function(x) log((1+sqrt(1-x^2))/x)-sqrt(1-x^2)

x0 = seq(1/2,1,by = .0001)

x1 = f(x0[1:(length(x0)-1)])

x2 = f(x0[2:(length(x0))])

x = x2 - x1

arc = sum(sqrt(.0001^2+x^2))

arc

###练习题2，求弧长。x = 1/4*y^2-1/2*log(y)  1<=y<=exp(1)

#使用循环结构求解

f = function(y) 1/4*y^2-1/2*log(y)

plot(f,xlim = c(1,exp(1)))

ub = exp(1)

lb = 1

m = 10000

step = (ub-lb)/m

y = 1

s = 0

for(i in 1:m){

s = s + sqrt(step^2 + (f(y+step)-f(y))^2)

y = y+step

}

cat(s)

#使用矩阵方法求解

f = function(y) 1/4*y^2-1/2*log(y)

y0 = seq(1,exp(1),by = .0001)

y1 = f(y0[1:(length(y0)-1)])

y2 = f(y0[2:(length(y0))])

y = y2 - y1

arc = sum(sqrt(.0001^2+y^2))

arc

 

###定积分在几何上的应用——求面积

f1 = function(x) sqrt(8-x^2)

f2 = function(x) .5*x^2

lb = 0

ub = 2

m = 10000

s = 0

step = (ub-lb)/m

x = 0

for(i in 1:m){

s = s + (f1(x)-f2(x))*step

x = x+step

}

s = 2*s

s

s0 = 8*pi-s

s0

 

 

###定积分在几何上的应用——求体积

f1 = function(x) x^2

f2 = function(x) sqrt(x)

lb = 0

ub = 1

m = 10000

v = 0

step = (ub-lb)/m

x = 0

for(i in 1:m){

sd = pi*(f2(x))^2-pi*(f1(x))^2

v = v+ sd*step

x = x+step

}

v