c语言和matlab语言实现牛顿迭代法解非线性方程和非线性方程组

#include<stdio.h>
#include<stdlib.h>  
#include<math.h>  
#include<float.h>
#include<time.h>

#define PI 3.14159265358979323846  /* pi */
#define ε 1.0e-12
int main()
{
    double x0 = PI;//取的初始值
    double x1 = 0.0;//有x0算出的x1，初始值先给定0
    double fx = 0.0;//f（x）
    double fxp = 0.0;//f（x）的导数
    double faix = 0.0;//计算结果,牛顿迭代格式 faix =x-(fx/fxp)
    int i = 0;//迭代计算次数
    while (fabs((x0 - x1) / x1) > ε)//判断两次迭代变量之间相对误差与精度比较
    {
        x1 = x0;//用x1存放上次的x0
        fx = sin(x0) - x0 / 2;
        fxp = cos(x0) - 0.5;
        faix = x0 - fx / fxp;
        x0 = faix;//将迭代后的结果赋给上次代入值变量，供下一次代入使用
        i++;//计算次数
        printf("第%d次迭代，迭代结果为：,%+.12e  
", i, x1);
    }
}

题目：计算sinx=x/2的根。

分析：newton法在大范围的收敛定理：

函数f（x）在区间[a，b]上存在二阶连续导数，且满足4个条件：

1.　　f（a）*f（b)<0

2.　　当x属于[a,b]时，函数的导数值不等于零。

3.　　当x属于[a,b]时，函数的二阶导数值保号。

4.　　a-f（a）/f'（a）<=b,且b-f（b）/f'（b）<=a

 计算结果：

---

matlab求解非线性方程：

，x=[pi/2,pi] 。

clc;
clear all;
close all;
%% 绘图
ezplot('sin(x)-x/2')
hold on;
ezplot('sin(x)')
hold on;
ezplot('x/2')
hold on;
ezplot('y=0*x')
legend('f(x)=sin(x)-x/2','sin(x)','x/2')
title('求解非线性方程')
%% 牛顿迭代法
fx=@(x)sin(x)-x/2.0;%定义 f(x)=sin(x)-x/2 匿名函数
DfxDx=@(x) cos(x)-1/2.0;% 定义f'(x)
epsilonT=1e-12;%收敛判断标准：相对误差
x0=pi;%给初值
n=0;%迭代次数
while 1
    x1=x0-fx(x0)/DfxDx(x0);
    epsilon=abs((x1-x0)/x1);
    x0=x1;
    n=n+1;
    disp(['第' num2str(n) '次迭代，x=', num2str(x1,15)]);
    if epsilon<epsilonT
        break;%跳出while循环
    end
end
disp('end')

第1次迭代，x=2.0943951023932
第2次迭代，x=1.91322295498104
第3次迭代，x=1.89567175194481
第4次迭代，x=1.89549428525543
第5次迭代，x=1.89549426703398
第6次迭代，x=1.89549426703398
end

---

 1stopt解非线性方程：

，x=[pi/2,pi] 。

NewCodeBlock"求解非线性方程";
//
AlgorithmOption =[1,0,1.00E-30,1000,20,30,50,20,0,0.15,1];
Parameter x=[pi/2,pi];
Function  sin(x)=x/2;

求解非线性方程

====== 结果 ======

迭代数: 21
计算用时(时:分:秒:微秒): 00:00:00:119
计算结束原因: 达到收敛判断标准
优化算法: 通用全局优化算法(UGO1)
函数表达式 sin(x)-(x/2) = 0
目标函数值(最小): 0
x: 1.89549426703398

====== 计算结束 ======

---

 准牛顿方法解非线性方程：sin(x)=x/2,x=[pi/2,pi] 

https://zhuanlan.zhihu.com/p/101077902

%% qusi-newton 准牛顿（割线法，不用求导数，用割线斜率代替切线）
clc;
clear all;
close all;
f=@(x)sin(x)-x/2.0;%定义 f(x)=sin(x)-x/2 匿名函数
epsilonT=1e-12;%收敛判断标准：相对误差
x0=pi/2;%给初值
x1=pi/2+0.1;
n=0;%迭代次数
while 1
%   x2=x0-f(x0)*(x1-x0)/(f(x1)-f(x0));
     x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0));%等同于上句，都行
    epsilon=abs((x1-x0)/x1);
    x0=x1;
    x1=x2;
    n=n+1;
    disp(['第' num2str(n) '次迭代，x=', num2str(x1,15)]);
    if epsilon<epsilonT
        break;%跳出while循环
    end
end
disp('qusi-newton end')

第1次迭代，x=1.96101103612234
第2次迭代，x=1.88595391153747
第3次迭代，x=1.89514618938836
第4次迭代，x=1.89549620158427
第5次迭代，x=1.89549426664428
第6次迭代，x=1.89549426703398
第7次迭代，x=1.89549426703398
第8次迭代，x=1.89549426703398
qusi-newton end

---

非线性多元方程组用牛顿法求解：

% https://zhuanlan.zhihu.com/p/146371408
% https://zhuanlan.zhihu.com/p/63103354 知乎代码
clc;
clear all;
close all;
x0=[1 2];
eps=1e-12;
[allx,ally,r,n]=mulNewton(fun,x0,eps);

disp(['迭代' num2str(n) '次，'  'x1=' num2str(eval(r(1)),100)  ',x2=' num2str(eval(r(2)),100) ])% 0.19808577588668504464406945200284
% disp(r)
disp('newton end')
%% 子函数区域
function [allx,ally,r,n]=mulNewton(F,x0,eps)
  if nargin==2
    eps=1.0e-4;
  end
  x0 = transpose(x0);
  Fx = subs(F,transpose(symvar(F)),x0);%将初始点代入方程组
  var = transpose(symvar(F));
  dF = jacobian(F,var);%求雅克比矩阵
  dFx = subs(dF,transpose(symvar(F)),x0);%将当前解带入雅克比矩阵
  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=2;
    for i=1:k
      x(i)=sym (['x',num2str(i)]);
    end 

  f(1)=0.5*sin(x(1))+0.1*cos(x(1)*x(2))-x(1);
  f(2)=0.5*cos(x(1))-0.1*cos(x(2))-x(2);
end

---

1stopt求解非线性方程组：

NewCodeBlock"求解非线性方程组";
//https://zhuanlan.zhihu.com/p/63103354
//https://zhuanlan.zhihu.com/p/146371408
AlgorithmOption =[1,0,1.00E-30,1000,20,30,50,20,0,0.15,1];
Parameter x(2);
Function  0.5*sin(x1)+0.1*cos(x1*x2)-x1=0;
          0.5*cos(x1)-0.1*cos(x2)-x2=0;

求解非线性方程组

====== 结果 ======

迭代数: 21
计算用时(时:分:秒:微秒): 00:00:00:429
计算结束原因: 达到收敛判断标准
优化算法: 通用全局优化算法(UGO1)
函数表达式 1: 0.5*sin(x1)+0.1*cos(x1*x2)-x1-0 = 0
2: 0.5*cos(x1)-0.1*cos(x2)-x2-0 = 0
目标函数值(最小): 0
x1: 0.198085775886685
x2: 0.398040303134032

====== 计算结束 ======

---

多元线性方程组求解：https://zhuanlan.zhihu.com/p/128577559

高斯列主元消去法：

% 高斯列主元消去法  https://zhuanlan.zhihu.com/p/128577559
% 先选择每列元素绝对值最大值放通过行变换放在主对角线上，之后将矩阵A化为上三角矩阵，然后回代求解线性方程组Ax=b。
clc;
clear all;
close all;
% A=[2 8 2;
%     1 6 -1;
%     2 1 2];%https://wenku.baidu.com/view/99bf58cf050876323112120c.html
% b=[14 13 5]';
A=[0.001 2 3;
    -1 3.712 4.623;
    -2 1.072 5.643];%https://wenku.baidu.com/view/99bf58cf050876323112120c.html
b=[1
    2 
    3 ];
x = gaussPivotScale(A, b)
disp('end');

function x = gaussPivotScale(A, b)

%GAUSSPIVOTSCALE  It uses Gauss elimination with partial pivoting
% and row scaling to solve the linear system Ax = b.
%   A : It is the n-by-n coefficient matrix.
%   b : It is the n-by-1 result vector.
%   x : It is the n-by-1 solution vector.

tic
format long

% judge wether there is a solution or not from the linear system
r_A = rank(A);
r_Ab = rank([A, b]);

if r_A ~= r_Ab
    disp('This linear system is no solution');
end

% Elimination

[row, ~] = size(A);
C = [A, b];

for k = 1 : row-1

% style of table
fprintf('This is%2dth transformation 
 [A | b] = 
', k);

    M = max(abs(C(k:end, k:end-1)), [], 2);  % find maximum in kth row,A所有行每行元素绝对值最大值
    a = abs(C(k:end, k));  % each value is absoluted value in kth column，所有行第k列绝对值
    [~, id] = max(a ./ M); % find row with maximum
    
    if id > k
        C([k, id], :) = C([id, k], :);  % pivot rows
    end
 
    mult = C(k+1:row, k) / C(k, k);  % construct multipliers
    
    % creat mesh, and improve speed
    [C_k, mult] = meshgrid(C(k, :), mult);  
    C(k+1:row, :) = C(k+1:row, :) - C_k .* mult;
    
disp(C)
end

% Back Substitution
[row, ~] = size(C); 

for k = row : -1 : 1
    C(k, :) = C(k, :) ./ C(k, k);
    C(1:k-1, row+1) = C(1:k-1, row+1) - C(1:k-1, k) * C(k, row+1);
end

x = C(:, end);

toc

end

---