數值計算——求解非線性方程組

        上一篇是求解的非線性方程，但是并不能簡單的擴充到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