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

我就是我 2022-06-16 00:41 294阅读 0赞

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

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

## 1、牛顿法 ##

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

![20170430223248681][]

算法：

![SouthEast][]

代码：

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、割线法 ##

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

![SouthEast 1][]

代码：

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

[20170430223248681]: /images/20220616/8b2f49b831884f68949922a3490af907.png
[SouthEast]: /images/20220616/a5c021df5d8247eb8bb970210eac7842.png
[SouthEast 1]: /images/20220616/a1fa58b4919f4998af55f9c27719be4e.png