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  \n", 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 \n [A | b] = \n', 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

 

 

 

---