R语言编程数学分析重读微积分微分学原理运用

x = seq(0,0.1,0.01)
y = runif(11,0,1)
plot(x,y)
lines(x,y)

x = seq(0,0.01,0.00001)
y = runif(1001,0,1)
plot(x,y)

x = seq(-10,10)
y = abs(x)
plot(x,y,type='l')

> x = seq(-1,1,0.1)
> y = sin(x)
> y1 = cos(x)
> xEnd = x+0.1
> yEnd = y+y1*0.1
> plot(x,y)
> for(i in seq_along(x)){
+ lines(c(x[i],xEnd[i]),c(y[i],yEnd[i]),col="red")
+ }

x = seq(-5,5,0.1)
y = cos(x)
x1 = seq(-5,5,0.1)
end = length(x1)
y1 = (sin(x1[2:end])-sin(x1[1:end-1]))/0.1
x5 = seq(-5,5,0.5)
end = length(x5)
y5 = (sin(x5[2:end])-sin(x5[1:end-1]))/0.5
x10 = seq(-5,5,1)
end = length(x10)
y10 = (sin(x10[2:end])-sin(x10[1:end-1]))/0.5
plot(x,y,type="l",col="red")
lines(x1[2:length(x1)],y1)
lines(x5[2:length(x5)],y5)
lines(x10[2:length(x10)],y10)

# 差商算法,x,y为同等长度的数组
diffDiv<-function(x,y){
    end = length(x)
    ind = 2:end #索引
    return(
        if(end==1) y[1]
        else (diffDiv(x[ind],y[ind])
            -diffDiv(x[ind-1],y[ind-1]))/(x[end]-x[1])
    )
}

kDiffDiv <- function(x,y,k){
    len <- length(x)
    if(len<k) return(0)
    d<-rep(0,len-k)
    for(i in 1:(len-k))
        d[i] <- diffDiv(x[i:(i+k)],y[i:(i+k)])
    return(d)
}

> plot(x,y)
> k = 1
> lines(x[1:(end-k)],kDiffDiv(x,y,k),col="red")
> k = 2
> lines(x[1:(end-k)],kDiffDiv(x,y,k),col="green")
> k = 3
> lines(x[1:(end-k)],kDiffDiv(x,y,k),col="blue")
> k = 4
> lines(x[1:(end-k)],kDiffDiv(x,y,k),col="purple")
> k = 5
> lines(x[1:(end-k)],kDiffDiv(x,y,k),col="yellow")
> k = 6
> lines(x[1:(end-k)],kDiffDiv(x,y,k),col="gray")

polyMulti<-function(x){
    n = length(x)
    if(n<2) return(c(-x,1))
    omega = rep(0,n+1)
    omega[1] = - x[1]
    omega[2] = 1
    for(i in 2:n){
        omega[2:i] = -x[i]*omega[2:i]*+omega[1:(i-1)]
        omega[1] = -x[i]*omega[1]
        omega[i+1] = 1
    }
    return(omega)
}
Newton<-function(x,y){
    n = length(x)
    N = rep(0,n+1)
    N[1]=y[1]
    for(i in 1:n){
        omega = polyMulti(x[1:i])
        d = kDiffDiv(x[1:i],y[1:i],i-1)
        N[1:i] = N[1:i] + d*omega[1:i]
        N[i+1] = d*omega[i+1]
    }
    return(N)
}

x = sort(rnORM(10))
y = x^5+2*x^4
N = Newton(x,y)
Y = y*0
for(i in 1:length(Y))
    for(j in 1:length(N))
        Y[i] = Y[i]+N[j]*x[i]^(j-1)
plot(x,y)
lines(x,Y)

x = matrix(data=seq(-5,4.95,0.05),nrow=200,ncol=200)
y = t(x)
z = 1-(x^2+y^2)
library(rgl)
persp3D(x,y,z,col="red")
N = 1000
x0 = rnorm(N)
y0 = rnorm(N)
z0 = 1-(x0^2+y0^2)
x1 = x0+3
z1 = -2*x0*3+z0
for(i in 1:N)
    lines3d(c(x0[i],x1[i]),c(y0[i],y0[i]),c(z0[i],z1[i]),col="green")

x = matrix(data=seq(-5,4.95,0.05),nrow=200,ncol=200)
y = t(x)
z = -(x^2+y^2)
library(rgl)
persp3d(x,y,z,col="red")
theta = seq(-pi,pi,0.01)
x0 = cos(theta)
y0 = sin(theta)
z0 = 1-(x0^2+y0^2)
x1 = x0*0
y1 = y0*0
z1 = z0+2
for(i in 1:N)
    lines3d(c(x0[i],x1[i]),c(y0[i],y1[i]),c(z0[i],z1[i]),col="green")

theta = seq(-pi,pi,0.1)
x = cos(theta)
y = sin(theta)
plot(x,y,type='l',col='red')
x1 = x*0
y1 = (x1-x)/(-2*x)*(-2*y)+y
for(i in 1:length(theta))
    lines(c(x[i],x1[i]),c(y[i],y1[i]),col='green')

library(rgl)
N = 1500
x<-rnorm(N)
y<-rnorm(N)
z<-1-x^2-y^2
for(i in 1:N){
    lines3d(c(x[i],3*x[i]),c(y[i],3*y[i]),c(z[i],z[i]+1),col='red')
    if(y[i]>0.1)
        lines3d(c(x[i],x[i]),c(y[i],y[i]-1/y[i]/2),c(z[i],z[i]+1),col='green')
    if(x[i]>0.1)
    lines3d(c(x[i],x[i]-1/x[i]/2),c(y[i],y[i]),c(z[i],z[i]+1),col='green')
    lines3d(c(x[i],x[i]*(1-2*z[i])/(2-2*z[i])),c(y[i],y[i]*(1-2*z[i])/(2-2*z[i])),c(z[i],z[i]+1),col='green')
}

library(imager)
img = load.image("lena.jpg")
dim(img)
[1] 512 512   1   3
gray = grayscale(img)
par(mfrow=c(1,2))
plot(img)
plot(gray)

dim(gray)
[1] 512 512   1   1
mat = array(gray,dim=c(512,512))
mat_x = diff(mat,1)
mat_y = t(diff(t(mat),1))
par(mfrow=c(1,2))
image(mat_x)
image(mat_y)

a=mat[256,]
par(mfrow=c(1,2))
plot(a,type='l')
plot(diff(a,1),type='l')

theta = c(pi/3,pi/4,pi/5,pi/6)
par(mfrow=c(2,2))
for(i in 1:4)
image(mat_x[,1:511]*cos(theta[i])+mat_y[1:511,]*sin(theta[i]))