北邮Matlab选修 Week 5 牛顿迭代法

问题描述

用牛顿法求解二元方程组的根。

数学模型

在$(x_k,y_k)$处将将$f_1$和$f_2$分别进行泰勒级数展开，并取线性部分作为近似，如下：
令泰勒级数展开的近似式等于零，得到线性方程组：
表示为矩阵形式：
求解这个线性方程组，可以得到迭代增量，进而更新将$x$和$y$的值：
重复这个过程，直到满足收敛条件（例如，两次迭代的差的绝对值小于某个预定的容忍度）或达到最大迭代次数。

程序设计

 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


%newton_method.m
function X=newton_method()

X=[1;1];  %取初始解
N=100;    %最大迭代次数
e=1e-5;   %设置误差限

syms x y;
func1=x^2*cos(2*x)+y^2*sin(2*y)-1; %定义函数1
func2=x^3+y^3-6*cos(2*x*y)+1;      %定义函数2

d_func1_1=diff(func1,x);
d_func2_1=diff(func2,x);
d_func1_2=diff(func1,y);
d_func2_2=diff(func2,y);


f_out=f(X,func1,func2);

for i=1:N
    d_f_out=d_f(X,d_func1_1,d_func2_1,d_func1_2,d_func2_2); %计算雅各比矩阵
    det_x=-inv(d_f_out)*f_out;                              %计算x的变化量
    X=X+det_x;                                              %x更新
    f_out=f(X,func1,func2);                                 %计算非线性方程组
    if norm(det_x)<e
        break;
    end
end
X=double(X);
end
function [out] = f(X,func1,func2)
syms x y;
out(1,1)=subs(func1,{x,y},{X(1,1),X(2,1)});
out(2,1)=subs(func2,{x,y},{X(1,1),X(2,1)});
end
function d_out= d_f(X,d_func1_1,d_func2_1,d_func1_2,d_func2_2)
syms x y;
d_out(1,1)=double(subs(d_func1_1,{x,y},{X(1,1),X(2,1)}));
d_out(2,1)=double(subs(d_func2_1,{x,y},{X(1,1),X(2,1)}));
d_out(1,2)=double(subs(d_func1_2,{x,y},{X(1,1),X(2,1)}));
d_out(2,2)=double(subs(d_func2_2,{x,y},{X(1,1),X(2,1)}));
end
end

计算结果

1
2
3
4
5
6


>> newton_method

ans =

    0.6219
    0.9685

第五周作业-牛顿迭代法

BUPT大一下Matlab选修

北邮Matlab选修 Week 5 牛顿迭代法

问题描述

数学模型

程序设计

计算结果