简单迭代法、牛顿法、弦割法、布洛依登法求解方程或方程组【Matlab】

142 阅读 0 评论 94 点赞

我是靠谱客的博主怡然画笔，这篇文章主要介绍简单迭代法、牛顿法、弦割法、布洛依登法求解方程或方程组【Matlab】，现在分享给大家，希望可以做个参考。

利用迭代法求解定非线性方程及方程组，使得误差不超过10^(-8)。同时应用迭代加速技术，提交迭代运算效率。

此题需要用到的MATLAB代码及附录：

附录6 二分法作根的隔离

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
% 附录6 二分法作根的隔离
%% 二分法作根的隔离
clear                             %清除变量
clc                               %清除命令行窗口代码
format long
aa=input('n请输入自变量x的区间：n');
x=[aa(1):0.1:aa(2)];
y=fun(x);
plot(x,y);
hold on
ezplot('0');
xlim([aa(1),aa(2)]);
k=0 ;                           %二分法迭代次数
e=aa(2)-aa(1);                 %区间长度大小
a=aa(1);b=aa(2);c=(a+b)/2;
err1=1e-4; err2=1e-4;
while e>=err1 && abs(fun(c))>=err2
    if fun(a)*fun(c)<0
        b=c;
    end
    if fun(a)*fun(c)>0
        a=c;
    end
    if fun(a)*fun(c)==0
        a=c;b=c;
    end
    c=(a+b)/2;
    e=e/2;
    k=k+1;
end
xx=c,                             %二分法求得的近似解
fid=fopen('data_xx1.txt','w');
fprintf(fid,'%.4f',xx);         %保存二分法求得的近似解（满足1e-4精度要求）
fclose(fid);
fid=fopen('data_xx2.txt','w');
fprintf(fid,'%.4f',b);          %保存二分法求得的近似解（不满足精度要求）用于两点弦割法
fclose(fid);

%% 应创建一个“fun.m”函数文件与第3节运行程序处于同一文件夹
function y = fun(x)       %定义目标函数
y = x.^6-5*x.^5+3*x.^4+x.^3-7*x.^2+7*x-20;
end

附录7 利用简单迭代法求解方程

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
% 附录7 利用简单迭代法求解方程
%% 利用简单迭代法求解方程
clear       %清除变量
clc         %清除命令行窗口代码
format long
err1=input('n请输入最大允许误差：n');
k=0 ;                   %迭代次数
x(1)=load('data_xx1.txt');
e=1.0+err1;
while abs(e)>=err1
    k=k+1;
    x(k+1)=(5*x(k)^5-3*x(k)^4-x(k)^3+7*x(k)^2-7*x(k)+20)^(1/6);
    e=x(k+1)-x(k);
end
xx=x(k+1)                %简单迭代法求得的近似解
fid=fopen('data_xx_results.txt','w');
fprintf(fid,'%.8f',xx);
fclose(fid);

附录8 利用牛顿法求解方程

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
% 附录8 利用牛顿法求解方程
%% 利用牛顿法求解方程
clear         %清除变量
clc           %清除命令行窗口代码
format long
err1=input('n请输入最大允许误差：n');
x(1)=load('data_xx1.txt');
e=1.0+err1;
syms t
fun1 =@(t) t.^6-5*t.^5+3*t.^4+t.^3-7*t.^2+7*t-20;
fun2 = diff(fun1(t));
fun2 = matlabFunction(fun2);
max1=1000;   %最大迭代次数
y(1)=fun1(x(1));
for k = 1: max1
    x(k+1) = x(k)-fun1(x(k))/fun2(x(k));     %Newton法的迭代公式
    y(k+1) = fun1(x(k+1));
    e = abs(x(k+1)-x(k));      %误差区间
    if e<err1
        break
    end
end
xx=x(k+1)                %牛顿迭代法求得的近似解
figure;
plot([1:k+1],y)
xlabel('X');             %横轴为迭代次数
ylabel('Y');             %纵轴为方程的值f(x)
grid;

附录9 利用两点弦割法求解方程

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
% 附录9 利用两点弦割法求解方程
%% 利用两点弦割法求解方程
clear         %清除变量
clc           %清除命令行窗口代码
format long
err1 =input('n请输入最大允许误差：n');
x(1)=load('data_xx.txt');  y(1) = fun(x(1));
x(2)=x(1)+1e-5;    y(2) = fun(x(2));
max1=1000;   %最大迭代次数
for k = 2: max1
    x(k+1) = x(k)-fun(x(k))*(x(k)-x(k-1))/(fun(x(k))-fun(x(k-1)));     %两点弦割法的迭代公式
    y(k+1) = fun(x(k+1));
    e = abs(x(k+1)-x(k));      %误差区间
        if  e < err1
            break
        end
end
xx=x(k+1)                %牛顿迭代法求得的近似解
figure;
plot([1:k+1],y)
xlabel('迭代步数');
ylabel('误差');
grid;

附录10 利用单点弦割法求解方程

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
% 附录10 利用单点弦割法求解方程
%% 利用单点弦割法求解方程
clear         %清除变量
clc           %清除命令行窗口代码
format long
err1 =input('n请输入最大允许误差：n');
x(1)=load('data_xx.txt');  y(1) = fun(x(1));
x(2)=x(1)+1e-2;   y(2) = fun(x(2));
max1=1000;   %最大迭代次数
for k = 2: max1 
    x(k+1) = x(k)-fun(x(k))*(x(k)-x(1))/(fun(x(k))-fun(x(1)));     %两点弦割法的迭代公式
    y(k+1) = fun(x(k+1));
    e = abs(x(k+1)-x(k));      %误差区间
        if  e < err1
            break
        end
end
xx=x(k+1)                %牛顿迭代法求得的近似解
figure;
plot([1:k+1],y)
xlabel('迭代步数');
ylabel('误差');
grid;

附录11 利用牛顿法求解非线性方程组

复制代码

% 附录11 利用牛顿法求解非线性方程组
%% 利用牛顿法求解非线性方程组
clear         %清除变量
clc           %清除命令行窗口代码
format long
err1 = input('n请输入最大允许误差：n');
x0 = input('n请输入给定初始向量(元素之间逗号间隔)：n');
[allx,ally,xx,n]=NNewton(fun,x0,err1)
xx1=double(xx)
fid=fopen('data_xx.txt','w');
fprintf(fid,'%.8fn',xx1');    %将方程组数值解存于'data_xx.txt'文件夹中
fclose(fid);

function [allx,ally,r,n]=NNewton(F,x0,eps)
  if nargin==2              %函数调用中给定函数输入参数的数目为2
    eps=1.0e-8;
  end
  x0 = transpose(x0);       %transpose转置向量或矩阵
  Fx = subs(F,transpose(symvar(F)),x0);
  var = transpose(symvar(F));
  dF = jacobian(F,var);
  dFx = subs(dF,transpose(symvar(F)),x0);
  n=dFx;
  r=x0-inv(dFx)*Fx';
  n=1;tol=1;
  N=100;
  symx=length(x0);
  ally=zeros(symx,N);
  allx=zeros(symx,N);

while tol>eps
    x0=r;
    Fx = subs(F,transpose(symvar(F)),x0);
    dFx = subs(dF,transpose(symvar(F)),x0);
    r=vpa(x0-inv(dFx)*Fx');          %牛顿法迭代
    tol=norm(r-x0)
    if(n>N)
        disp('迭代步数太多，可能不收敛。');
        break;
    end
    allx(:,n)=x0;
    ally(:,n)=Fx;
    n=n+1;
  end
end

function f = fun(x)               %非线性方程组，根据题目或实际情况进行更改
  k=3;                              % k为所求解方程个数，同时也是自变量的个数
  for i=1:k                        %将自变量变成符号
      x(i)=sym (['x',num2str(i)]);
  end
  f(1)=x(1)^2+x(2)^2+x(3)^2-1;
  f(2)=2*x(1)^2+x(2)^2-4*x(3)^2;
  f(3)=3*x(1)^2-4*x(2)^2+x(3)^2;
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
% 附录11 利用牛顿法求解非线性方程组
%% 利用牛顿法求解非线性方程组
clear         %清除变量
clc           %清除命令行窗口代码
format long
err1 = input('n请输入最大允许误差：n');
x0 = input('n请输入给定初始向量(元素之间逗号间隔)：n');
[allx,ally,xx,n]=NNewton(fun,x0,err1)
xx1=double(xx)
fid=fopen('data_xx.txt','w');
fprintf(fid,'%.8fn',xx1');    %将方程组数值解存于'data_xx.txt'文件夹中
fclose(fid);

function [allx,ally,r,n]=NNewton(F,x0,eps)
  if nargin==2              %函数调用中给定函数输入参数的数目为2
    eps=1.0e-8;
  end
  x0 = transpose(x0);       %transpose转置向量或矩阵
  Fx = subs(F,transpose(symvar(F)),x0);
  var = transpose(symvar(F));
  dF = jacobian(F,var);
  dFx = subs(dF,transpose(symvar(F)),x0);
  n=dFx;
  r=x0-inv(dFx)*Fx';
  n=1;tol=1;
  N=100;
  symx=length(x0);
  ally=zeros(symx,N);
  allx=zeros(symx,N);

  while tol>eps
    x0=r;
    Fx = subs(F,transpose(symvar(F)),x0);
    dFx = subs(dF,transpose(symvar(F)),x0);
    r=vpa(x0-inv(dFx)*Fx');          %牛顿法迭代
    tol=norm(r-x0)
    if(n>N)
        disp('迭代步数太多，可能不收敛。');
        break;
    end
    allx(:,n)=x0;
    ally(:,n)=Fx;
    n=n+1;
  end
end

function f = fun(x)               %非线性方程组，根据题目或实际情况进行更改
  k=3;                              % k为所求解方程个数，同时也是自变量的个数
  for i=1:k                        %将自变量变成符号
      x(i)=sym (['x',num2str(i)]);
  end
  f(1)=x(1)^2+x(2)^2+x(3)^2-1;
  f(2)=2*x(1)^2+x(2)^2-4*x(3)^2;
  f(3)=3*x(1)^2-4*x(2)^2+x(3)^2;
end

附录12 利用弦割法求解非线性方程组

复制代码

% 附录12 利用弦割法求解非线性方程组
%% 利用弦割法求解非线性方程组(对应7.3第一个)
clear         %清除变量
clc           %清除命令行窗口代码
err1 = input('n请输入最大允许误差：n');                 % 求解精度
x0 = input('n请输入给定初始向量(元素之间逗号或空格间隔)：n');
x0=x0';
format long

f = @(x)[x(1).^2+x(2).^2+x(3).^2-1;    %待求函数，可以更改
    2*x(1).^2+x(2).^2-4*x(3);
3*x(1).^2-4*x(2).^2+x(3).^2];

h = 0.01;      % 步长
n=1;           %迭代次数
z = [f(x0 + h*[1; 0; 0]), f(x0 + h*[0; 1; 0]), f(x0 + h*[0; 0; 1])]f(x0);
x1 = x0 + h*z/(sum(z) - 1);
xx(:,n)=x1;
while norm(x1 - x0,2) >  err1       %弦割法迭代求解
    x0 =  x1;   n=n+1; 
    z = [f(x0 + h*[1; 0; 0]), f(x0 + h*[0; 1; 0]), f(x0 + h*[0; 0; 1])]f(x0);
    x1 = x0 + h*z/(sum(z) - 1);
    xx(:,n)=x1;
end
csvwrite('data_xx.txt', x1);

%% 利用弦割法求解非线性方程组(对应7.3第二个)
clear         %清除变量
clc           %清除命令行窗口代码
err1 = input('n请输入最大允许误差：n');                 % 求解精度
x0 = input('n请输入给定初始向量(元素之间逗号或空格间隔)：n');
x0=x0';
format long

f = @(x)[cos(x(1).^2+0.4*x(2))+x(1)^2+x(2)^2-1.6;    %待求函数，可以更改
1.5*x(1).^2-1/0.36*x(2).^2-1.0];

h = 0.01;      % 步长
n=1;           %迭代次数
z = [f(x0 + h*[1; 0]), f(x0 + h*[0; 1]) ]f(x0);
x1 = x0 + h*z/(sum(z) - 1);
xx(:,n)=x1;
while norm(x1-x0,2) >  err1       %弦割法迭代求解
    x0 = x1;   n = n+1; 
    z = [f(x0 + h*[1; 0 ]), f(x0 + h*[0; 1])]f(x0);
    x1 = x0 + h*z/(sum(z) - 1);
    xx(:,n)=x1;
end
csvwrite('data_xx.txt', x1);

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
% 附录12 利用弦割法求解非线性方程组
%% 利用弦割法求解非线性方程组(对应7.3第一个)
clear         %清除变量
clc           %清除命令行窗口代码
err1 = input('n请输入最大允许误差：n');                 % 求解精度
x0 = input('n请输入给定初始向量(元素之间逗号或空格间隔)：n');
x0=x0';
format long

f = @(x)[x(1).^2+x(2).^2+x(3).^2-1;    %待求函数，可以更改
    2*x(1).^2+x(2).^2-4*x(3);
3*x(1).^2-4*x(2).^2+x(3).^2];

h = 0.01;      % 步长
n=1;           %迭代次数
z = [f(x0 + h*[1; 0; 0]), f(x0 + h*[0; 1; 0]), f(x0 + h*[0; 0; 1])]f(x0);
x1 = x0 + h*z/(sum(z) - 1);
xx(:,n)=x1;
while norm(x1 - x0,2) >  err1       %弦割法迭代求解
    x0 =  x1;   n=n+1; 
    z = [f(x0 + h*[1; 0; 0]), f(x0 + h*[0; 1; 0]), f(x0 + h*[0; 0; 1])]f(x0);
    x1 = x0 + h*z/(sum(z) - 1);
    xx(:,n)=x1;
end
csvwrite('data_xx.txt', x1);

%% 利用弦割法求解非线性方程组(对应7.3第二个)
clear         %清除变量
clc           %清除命令行窗口代码
err1 = input('n请输入最大允许误差：n');                 % 求解精度
x0 = input('n请输入给定初始向量(元素之间逗号或空格间隔)：n');
x0=x0';
format long

f = @(x)[cos(x(1).^2+0.4*x(2))+x(1)^2+x(2)^2-1.6;    %待求函数，可以更改
1.5*x(1).^2-1/0.36*x(2).^2-1.0];

h = 0.01;      % 步长
n=1;           %迭代次数
z = [f(x0 + h*[1; 0]), f(x0 + h*[0; 1]) ]f(x0);
x1 = x0 + h*z/(sum(z) - 1);
xx(:,n)=x1;
while norm(x1-x0,2) >  err1       %弦割法迭代求解
    x0 = x1;   n = n+1; 
    z = [f(x0 + h*[1; 0 ]), f(x0 + h*[0; 1])]f(x0);
    x1 = x0 + h*z/(sum(z) - 1);
    xx(:,n)=x1;
end
csvwrite('data_xx.txt', x1);

附录13 利用布洛依登法求解非线性方程组

复制代码

% 附录13 利用布洛依登法求解非线性方程组
%% 布洛依登法求解非线性方程组
% 注释部分对应第一组方程，非注释代码对应第二组方程
clear         %清除变量
clc           %清除命令行窗口代码
format long
err1 = input('n请输入最大允许误差：n');
x0 = input('n请输入给定初始向量(元素之间逗号或空格间隔)：n');
x0=x0';
[x,n,data]=broyden(x0,err1);
disp('计算结果为')
x
disp('迭代次数为')
n
data(:,1)=x0;

%抽取 data 中第一个变量数据画出曲线
subplot(2,1,1)
plot(data(1,:)),
title('x1 在迭代中的变化')
xlim([1.0 5.0]);ylim([0.90 1.10]);
%抽取 data 中的第二个变量数据 画出其变化曲线
subplot(2,1,2)
plot(data(2,:)),title('x2 在迭代中的变化')
xlim([1.0 5.0]);ylim([0.40 0.50]);
% subplot(3,1,1)
% plot(data(1,:)),
% title('x1 在迭代中的变化')
% xlim([1.0 16.0]);ylim([0.00 1.00]);
% %抽取 data 中的第二个变量数据 画出其变化曲线
% subplot(3,1,2)
% plot(data(2,:)),title('x2 在迭代中的变化')
% xlim([1.0 16.0]);ylim([0.00 1.00]);
%抽取第三个变量数据
% subplot(3,1,3)
% plot(data(3,:)),title('x3 在迭代中的变化')
% xlim([1.0 16.0]);ylim([0.00 1.00]);

%% 布洛依登法计算非线性方程组
function [x,n,data]=broyden(x0,tol)
%输入 x0 为迭代初值，tol 为误差容限 如果缺省 默认为 10 的-10 次方
%data 用来存放计算的中间数据便于计算收敛情况分析
if nargin==1
tol=1e-8;
end
% H0=df2(x0);
H0=df3(x0);
H0=inv(H0);
% x1=x0-H0*f2(x0);
x1=x0-H0*f3(x0);
n=1;
%设置初始误差 使之可以进入循环
wucha=0.1;
%循环迭代
while (wucha>tol)&(n<20) &(n<500)
wucha=norm(x1-x0);
dx=x1-x0;
% y=f2(x1)-f2(x0);
y=f3(x1)-f3(x0);
fenzi=dx'*H0*y;                   %是标量
H1=H0+(dx-H0*y)*(dx)'*H0/fenzi;

temp_x0=x0;
x0=x1;
% x1=temp_x0-H1*f2(temp_x0); %x1 的更新
x1=temp_x0-H1*f3(temp_x0); %x1 的更新

%更新 H 矩阵
H=H1;
n=n+1;
%data 用来存放中间数据
data(:,n)=x1;
end
x=x1;
end

% 目标方程组(对应7.3第一个)，注意方程组输出为列向量函数
function F = f2(x0)
x=x0(1);
y=x0(2);
z=x0(3);
ff1=x^2+y^2+z^2-1.0;             %待求方程组，可更改
ff2=2*x^2+y^2-4*z;
ff3=3*x^2-4*y^2+z^2;
F=[ff1;ff2;ff3];
end

%Jacobi 矩阵,该方程输出为 n X n 维方阵(对应7.3第一个非线性方程组)
function f=df2(x0)
x=x0(1);
y=x0(2);
z=x0(3);
f=[2*x 2*y 2*z      %目标方程组Jacobi矩阵，可以更改（第一列对x求导，第二列对y求导，依次类推）
4*x 2*y -4
6*x -8*y 2*z];
end

% 目标方程组(对应7.3第二个)，注意方程组输出为列向量函数
function F = f2(x0)
x=x0(1);
y=x0(2);
ff1=cos(x^2+0.4*y)+x^2+y^2-1.6;             %待求方程组，可更改
ff2=1.5*x^2-1/0.36*y^2-1;
F=[ff1;ff2];
end

% Jacobi 矩阵,该方程输出为 2 X 2 维方阵(对应7.3第二个非线性方程组)
function f=df2(x0)
x=x0(1);
y=x0(2);
f=[-sin(x^2+0.4*y)*2*x+2*x -sin(x^2+0.4*y)*0.4+2*y    %目标方程组Jacobi矩阵，可以更改（第一列对x求导，第二列对y求导）
3*x -2/0.36*y]; 
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
% 附录13 利用布洛依登法求解非线性方程组
%% 布洛依登法求解非线性方程组
% 注释部分对应第一组方程，非注释代码对应第二组方程
clear         %清除变量
clc           %清除命令行窗口代码
format long
err1 = input('n请输入最大允许误差：n');
x0 = input('n请输入给定初始向量(元素之间逗号或空格间隔)：n');
x0=x0';
[x,n,data]=broyden(x0,err1);
disp('计算结果为')
x
disp('迭代次数为')
n
data(:,1)=x0;

%抽取 data 中第一个变量数据画出曲线
subplot(2,1,1)
plot(data(1,:)),
title('x1 在迭代中的变化')
xlim([1.0 5.0]);ylim([0.90 1.10]);
%抽取 data 中的第二个变量数据 画出其变化曲线
subplot(2,1,2)
plot(data(2,:)),title('x2 在迭代中的变化')
xlim([1.0 5.0]);ylim([0.40 0.50]);
% subplot(3,1,1)
% plot(data(1,:)),
% title('x1 在迭代中的变化')
% xlim([1.0 16.0]);ylim([0.00 1.00]);
% %抽取 data 中的第二个变量数据 画出其变化曲线
% subplot(3,1,2)
% plot(data(2,:)),title('x2 在迭代中的变化')
% xlim([1.0 16.0]);ylim([0.00 1.00]);
%抽取第三个变量数据
% subplot(3,1,3)
% plot(data(3,:)),title('x3 在迭代中的变化')
% xlim([1.0 16.0]);ylim([0.00 1.00]);

%% 布洛依登法计算非线性方程组
function [x,n,data]=broyden(x0,tol)
%输入 x0 为迭代初值，tol 为误差容限 如果缺省 默认为 10 的-10 次方
%data 用来存放计算的中间数据便于计算收敛情况分析
if nargin==1
tol=1e-8;
end
% H0=df2(x0);
H0=df3(x0);
H0=inv(H0);
% x1=x0-H0*f2(x0);
x1=x0-H0*f3(x0);
n=1;
%设置初始误差 使之可以进入循环
wucha=0.1;
%循环迭代
while (wucha>tol)&(n<20) &(n<500)
wucha=norm(x1-x0);
dx=x1-x0;
% y=f2(x1)-f2(x0);
y=f3(x1)-f3(x0);
fenzi=dx'*H0*y;                   %是标量
H1=H0+(dx-H0*y)*(dx)'*H0/fenzi;

temp_x0=x0;
x0=x1;
% x1=temp_x0-H1*f2(temp_x0); %x1 的更新
x1=temp_x0-H1*f3(temp_x0); %x1 的更新

%更新 H 矩阵
H=H1;
n=n+1;
%data 用来存放中间数据
data(:,n)=x1;
end
x=x1;
end

% 目标方程组(对应7.3第一个)，注意方程组输出为列向量函数
function F = f2(x0)
x=x0(1);
y=x0(2);
z=x0(3);
ff1=x^2+y^2+z^2-1.0;             %待求方程组，可更改
ff2=2*x^2+y^2-4*z;
ff3=3*x^2-4*y^2+z^2;
F=[ff1;ff2;ff3];
end

%Jacobi 矩阵,该方程输出为 n X n 维方阵(对应7.3第一个非线性方程组)
function f=df2(x0)
x=x0(1);
y=x0(2);
z=x0(3);
f=[2*x 2*y 2*z      %目标方程组Jacobi矩阵，可以更改（第一列对x求导，第二列对y求导，依次类推）
4*x 2*y -4
6*x -8*y 2*z];
end

% 目标方程组(对应7.3第二个)，注意方程组输出为列向量函数
function F = f2(x0)
x=x0(1);
y=x0(2);
ff1=cos(x^2+0.4*y)+x^2+y^2-1.6;             %待求方程组，可更改
ff2=1.5*x^2-1/0.36*y^2-1;
F=[ff1;ff2];
end

% Jacobi 矩阵,该方程输出为 2 X 2 维方阵(对应7.3第二个非线性方程组)
function f=df2(x0)
x=x0(1);
y=x0(2);
f=[-sin(x^2+0.4*y)*2*x+2*x -sin(x^2+0.4*y)*0.4+2*y    %目标方程组Jacobi矩阵，可以更改（第一列对x求导，第二列对y求导）
3*x -2/0.36*y]; 
end