数值计算——求解非线性方程组

数值计算——求解非线性方程组

        上一篇是求解的非线性方程，但是并不能简单的扩展到n维的情况下。下面是使用牛顿法和割线法求解非线性方程组：

1、牛顿法

      需要求方程组的雅克比矩阵，即对每个变量求偏导：



      算法：



       代码：

package com.kexin.lab7;
import com.kexin.lab3.*;
/**
* 牛顿法求解非线性方程组
* @author KeXin
*
*/
public class Newton {
static final double tol = 0.00001;
/**
* 返回函数值的取反,这里处理的函数是：
* (x1+3)(x2^3-7)+18=0
* sin(x2*e^x1-1)=0
* @param x
* @return
*/
public static double[] Function(double[] x){
double[]result = new double[2];
double x1 = x[0];
double x2 = x[1];
//书上例题测试
//		result[0] = -(x1+2*x2-2);
//		result[1] = -(Math.pow(x1, 2)+4*Math.pow(x2, 2)-4);
//实验题目
result[0] = -((x1+3)*(Math.pow(x2, 3)-7)+18);
result[1] = -(Math.sin((x2*Math.pow(Math.E, x1)-1)));
return result;
}
/**
* 返回上面函数对应的Jacobi矩阵
* @param x
* @return
*/
public static double[][] FunctionJ(double[] x){
double[][] result = new double[2][2];
double x1 = x[0];
double x2 = x[1];
//书上例题测试
//		result[0][0] = 1;
//		result[0][1] = 2;
//		result[1][0] = 2*x1;
//		result[1][1] = 8*x2;
//实验题目
result[0][0] = Math.pow(x2, 3)-7;
result[0][1] = (3*x1+9)*Math.pow(x2, 2);
double e = Math.E;
double temp = Math.cos((x2*Math.pow(e, x1)-1))*Math.pow(e, x1);
result[1][0] = x2*temp;
result[1][1] = temp;
return result;
}
public static double[] NewtonF(double[]x0){
double[] f ;
double[][] j ;
double[] s;
double r;
double[] x = {0,1};
int time = 0;
while(true){
time++;
f = Function(x0);
if(Math.abs(f[0])<tol&&Math.abs(f[1])<tol){
break;
}
j = FunctionJ(x0);
DieDai.PrintArray("x", x0);
DieDai.PrintArray("-F", f);
DieDai.Print2DArray("J", j);
s = DieDai.SOR(j, f,0.2);
DieDai.PrintArray("s", s);
r = Math.sqrt((Math.pow(x[0]-x0[0], 2)+Math.pow(x[1]-x0[1],2)));
System.out.println("r:"+r);
x0 = AddFunction(x0, s);
}
System.out.println("time = "+time);
return x0;
}
public static void main(String[] args) {
//定义初值x0
//		double[] x0 = {1,2};	//书上例题测试
double[] x0 = {-0.5,1.4};	//实验题目
double[] result = NewtonF(x0);
DieDai.PrintArray("result", result);
DieDai.PrintArray("-F", Function(result));
}
/**
* 转置行矩阵
*
* @param a
* @return
*/
public static double[][] Transposition1(double[] a) {
int len = a.length;
double[][] result = new double[len][1];
for (int i = 0; i < len; i++) {
result[i][0] = a[i];
}

return result;
}
/**
* 一维数组矩阵加法
* @param x1
* @param x2
* @return
*/
public static double[] AddFunction(double[]x1,double[]x2){
int row = x1.length;
double[] result = new double[row];
for(int i = 0;i<row;i++){
result[i] = x1[i]+x2[i];
}
return result;
}

}

2、割线法

          割线法是牛顿法的变形之一，牛顿法是使用导数即割线逼近正确解，而割线法则是选择两点的连线即曲线的割线逼近正确解：



       代码：
broyden.m

function [time,x,f] = broyden()
x = [-0.5;1.4]
B=[1,2;2,6];
i = 0
time = 0
while i<1
time = time + 1
f=caculatef(x);
s=inv(B)*f;
s=-s;
x=x+s;
y=caculatef(x)-f;
B=B+inv((s'*s))*(y-B*s)*s';
if s(1,1)^2+s(2,1)^2<0.00001
i=2
end
end
end

Caculatef.m:

function [f] = caculatef(x)
f(1,1)=(x(1,1)+3)*(x(2,1)^3-7)+18;
f(2,1)=sin(x(2,1)*exp(x(1,1))-1);
end