最近接触到一些工学上的数学计算,琢磨着如何用R来轻松秒杀,总结出了一些方法。
1. 使用逆行列矩阵
举个例子,解下面的二元一次方程
2x+3y=8
x-2y=5
这个用最简单的线性代数就行,Av=b。
R的话用到solve
来求逆行列。
> A=matrix(c(2,1,3,-2),2,2)
> A
[,1] [,2]
[1,] 2 3
[2,] 1 -2
> b=c(8,5)
> b
[1] 8 5
> solve(A,b)
[1] 4.4285714 -0.2857143
三元方程的话也是一样的,只不过A多了一列而已。
2. 用nleqslv函数
R还有专门的求解方程的包nleqslv。当然不只是解方程,求最优值什么的也可以用这个包。安装完了包来解决同样的问题。
g <- function(z)
{ x = z[1]
y = z[2]
f1 = 2*x+3*y-8
f2 = x-2*y-5
c(f1,f2)
}
nleqslv(c(1, 1), g) # 制定起始值
可以求得同样的结果
> result<-nleqslv(c(1, 1), g)
> result$x
[1] 4.4285714 -0.2857143
当然这只是最简单的例子,这个包主要用于计算高次元非线形方程。
比方说
3x+xy=15
x*x+2y=20
这个时候解肯定不止一个,用R的话可以先把这两个函数可视化一下,然后用可交互界面鼠标点击两个函数的交点,给出初始值,求解。
x=seq(-10, 20, 0.5)
y=seq(-10,20,0.5)
f1=function(x,y) 3*x+x*y - 15
f2=function(x,y) x^2+y*2-20
z1 = outer(x, y, f1)
z2 = outer(x, y, f2)
contour(x, y, z1, levels=0,ylim=c(-10,20))
contour(x, y, z2, levels=0,col=2,add=T)
abline(h=0,col=8); abline(v=0,col=8)
grid()
#定义联列方程式
library(nleqslv)
g= function(z) {
x = z[1]
y = z[2]
f1 = 3*x+ x*y - 15 # 方程式1
f2 = x^2+y*2-20 # 方程式2
c(f1, f2)
}
#求解
p=locator(1); points(p,col=4) # 指定初始值
x0=p$x ; y0=p$y # 鼠标点击指定初始值
ans = nleqslv(c(x0,y0), g) #求解
ans$x
points(ans$x[1],ans$x[2],pch=20,col=2)
q=locator(1); # 用鼠标点击指定text的位置
text(q$x,q$y,paste("x=",round(ans$x[1],3),"y=",round(ans$x[2],3)))
以上