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

        上一篇是求解的非线性方程，但是并不能简单的扩展到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 = 0time = 0while 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    endendend

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