Newton法说明:
Newton法R程序:
例
{x21+x22−5=0(x1+1)∗x2−(3∗x1+1)=0
{
x
1
2
+
x
2
2
−
5
=
0
(
x
1
+
1
)
∗
x
2
−
(
3
∗
x
1
+
1
)
=
0
R程序
Newtons<-function(fun,x,ep=1e-5,it_max=100){#fun是函数,x表示初值,ep表示精度,itmax表示最大迭代次数
index<-0;k<-1
while(k<=it_max){
x1<-x;obj<-fun(x);
x<-x-solve(obj$j,obj$f)
norm<-sqrt((x-x1)%*%(x-x1))
if(norm<ep){
index<-1;break
}
k=k+1
}
obj<-fun(x)
list(root=x,it=k,index=index,Fval=obj$f)
}
#例x[1]^2+x[2]^2-5=0;(x[1]+1)*x[2]-(3*x[1]+1)=0
fun<-function(x){
f<-c(x[1]^2+x[2]^2-5,(x[1]+1)*x[2]-(3*x[1]+1))
j<-matrix(c(2*x[1],2*x[2],x[2]-3,x[1]+1),nrow = 2,byrow = T)
list(f=f,j=j)#f是函数值,j是雅可比行列式
}
Newtons(fun,c(0,1))
其中部分说明
- list对应$使用;function(函数)返回值就是函数的输出项(没有对应要赋值)
> a<-list(f=1,j=2)
> a$f
[1] 1
> a<-function(x){list(f=1,j=2)}
> a(1)$j
[1] 2
%*%表示转置相乘更多关于运算符详见
solve()函数
例求
⎧⎩⎨⎪⎪2x1+0x2+2x3=12x1+1x2+2x3=22x1+1x2+0x3=3 { 2 x 1 + 0 x 2 + 2 x 3 = 1 2 x 1 + 1 x 2 + 2 x 3 = 2 2 x 1 + 1 x 2 + 0 x 3 = 3
> A
[,1] [,2] [,3]
[1,] 2 0 2
[2,] 2 1 2
[3,] 2 1 0
> b=1:3
> b
[1] 1 2 3
> solve(A,b)
[1] 1.0 1.0 -0.5
即
A−1∗b
A
−
1
∗
b
也就是Newton法中要求的
J−1∗f(xk)
J
−
1
∗
f
(
x
k
)
具体 详见
- matrix()函数
把向量变成矩阵的函数
nrow
the desired number of rows.(行数)
byrow
logical. If FALSE (the default) the matrix is filled by columns, otherwise the matrix is filled by rows.(默认为F按照列排,T表示按照行排)
如:
> matrix(c(1,2,3,4),nrow = 2)
[,1] [,2]
[1,] 1 3
[2,] 2 4
> matrix(c(1,2,3,4),nrow = 2,byrow = T)
[,1] [,2]
[1,] 1 2
[2,] 3 4
练习题1
编写一个用newton法求解非线性方程根的函数,并求函数
x3−x−1=0
x
3
−
x
−
1
=
0
终止条件为
|xk−xk−1|<10−5
|
x
k
−
x
k
−
1
|
<
10
−
5
代码
#例x^3-x-1=0
fun<-function(x){
f<-x^3-x-1
j<-3*x
list(f=f,j=j)
}
x=1;#随便定个初值
it_max=100;index=0;k=1
while(k<=it_max){
x1=x;
x<-x-solve(fun(x)$j,fun(x)$f)
if(sqrt((x1-x)%*%(x1-x))<1e-5) {index=1;break}
k=k+1
}
list(root=x,it=k,index=index,fval=fun(x)$f)
#或者直接调用之前的Newtons函数
Newtons(fun,1)
结果
> #例x^3-x-1=0
> fun<-function(x){
+ f<-x^3-x-1
+ j<-3*x
+ list(f=f,j=j)
+ }
> x=1;#随便定个初值
> it_max=100;index=0;k=1
> while(k<=it_max){
+ x1=x;
+ x<-x-solve(fun(x)$j,fun(x)$f)
+ if(sqrt((x1-x)%*%(x1-x))<1e-5) {index=1;break}
+ k=k+1
+ }
> list(root=x,it=k,index=index,fval=fun(x)$f)
$`root`
[1] 1.324718
$it
[1] 5
$index
[1] 1
$fval
[1] 1.070586e-06
> Newtons(fun,1)
$`root`
[1] 1.324718
$it
[1] 5
$index
[1] 1
$Fval
[1] 1.070586e-06
练习题2
用Newton法求非线性方程组
{x21+x22−1=0x31−x2=0
{
x
1
2
+
x
2
2
−
1
=
0
x
1
3
−
x
2
=
0
取初始值
x(0)=(−0.8,0.6)T
x
(
0
)
=
(
−
0.8
,
0.6
)
T
,精度要求
ξ=10−3
ξ
=
10
−
3
R代码
#例 x_1^2+x_2^2-1=0\\x_1^3-x_2=0 精度10^(-3)
fun<-function(x){
f<- c(x[1]^2+x[2]^2-1,x[1]^3-x[2])
j<-matrix(c(2*x[1],2*x[2],3*x[1]^2,-1),nrow = 2,byrow = T)
list(f=f,j=j)
}
#直接调用之前的Newtons函数(因为精度变了,稍微改下)
Newtons<-function(fun,x,ep,it_max=100){#itmax表示最大迭代次数
index<-0;k<-1
while(k<=it_max){
x1<-x;obj<-fun(x);
x<-x-solve(obj$j,obj$f)
norm<-sqrt((x-x1)%*%(x-x1))
if(norm<ep){
index<-1;break
}
k=k+1
}
obj<-fun(x)
list(root=x,it=k,index=index,Fval=obj$f)
}
Newtons(fun,c(-0.8,0.6),1e-3)
结果
> #例 x_1^2+x_2^2-1=0\\x_1^3-x_2=0 精度10^(-3)
> fun<-function(x){
+ f<- c(x[1]^2+x[2]^2-1,x[1]^3-x[2])
+ j<-matrix(c(2*x[1],2*x[2],3*x[1]^2,-1),nrow = 2,byrow = T)
+ list(f=f,j=j)
+ }
> #直接调用之前的Newtons函数(因为精度变了,稍微改下)
> Newtons<-function(fun,x,ep,it_max=100){#itmax表示最大迭代次数
+ index<-0;k<-1
+ while(k<=it_max){
+ x1<-x;obj<-fun(x);
+ x<-x-solve(obj$j,obj$f)
+ norm<-sqrt((x-x1)%*%(x-x1))
+ if(norm<ep){
+ index<-1;break
+ }
+ k=k+1
+ }
+ obj<-fun(x)
+ list(root=x,it=k,index=index,Fval=obj$f)
+ }
> Newtons(fun,c(-0.8,0.6),1e-3)
$`root`
[1] 0.8260317 0.5636240
$it
[1] 6
$index
[1] 1
$Fval
[1] 4.106884e-07 9.268858e-07