c语言和matlab语言实现牛顿迭代法解非线性方程和非线性方程组
1 #include<stdio.h> 2 #include<stdlib.h> 3 #include<math.h> 4 #include<float.h> 5 #include<time.h> 6 7 #define PI 3.14159265358979323846 /* pi */ 8 #define ε 1.0e-12 9 int main() 10 { 11 double x0 = PI;//取的初始值 12 double x1 = 0.0;//有x0算出的x1，初始值先给定0 13 double fx = 0.0;//f（x） 14 double fxp = 0.0;//f（x）的导数 15 double faix = 0.0;//计算结果,牛顿迭代格式 faix =x-(fx/fxp) 16 int i = 0;//迭代计算次数 17 while (fabs((x0 - x1) / x1) > ε)//判断两次迭代变量之间相对误差与精度比较 18 { 19 x1 = x0;//用x1存放上次的x0 20 fx = sin(x0) - x0 / 2; 21 fxp = cos(x0) - 0.5; 22 faix = x0 - fx / fxp; 23 x0 = faix;//将迭代后的结果赋给上次代入值变量，供下一次代入使用 24 i++;//计算次数 25 printf("第%d次迭代，迭代结果为：,%+.12e \n", i, x1); 26 } 27 }
题目：计算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] 。
1 clc; 2 clear all; 3 close all; 4 %% 绘图 5 ezplot(\'sin(x)-x/2\') 6 hold on; 7 ezplot(\'sin(x)\') 8 hold on; 9 ezplot(\'x/2\') 10 hold on; 11 ezplot(\'y=0*x\') 12 legend(\'f(x)=sin(x)-x/2\',\'sin(x)\',\'x/2\') 13 title(\'求解非线性方程\') 14 %% 牛顿迭代法 15 fx=@(x)sin(x)-x/2.0;%定义 f(x)=sin(x)-x/2 匿名函数 16 DfxDx=@(x) cos(x)-1/2.0;% 定义f\'(x) 17 epsilonT=1e-12;%收敛判断标准：相对误差 18 x0=pi;%给初值 19 n=0;%迭代次数 20 while 1 21 x1=x0-fx(x0)/DfxDx(x0); 22 epsilon=abs((x1-x0)/x1); 23 x0=x1; 24 n=n+1; 25 disp([\'第\' num2str(n) \'次迭代，x=\', num2str(x1,15)]); 26 if epsilon<epsilonT 27 break;%跳出while循环 28 end 29 end 30 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] 。
1 NewCodeBlock"求解非线性方程"; 2 // 3 AlgorithmOption =[1,0,1.00E-30,1000,20,30,50,20,0,0.15,1]; 4 Parameter x=[pi/2,pi]; 5 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
1 %% qusi-newton 准牛顿（割线法，不用求导数，用割线斜率代替切线） 2 clc; 3 clear all; 4 close all; 5 f=@(x)sin(x)-x/2.0;%定义 f(x)=sin(x)-x/2 匿名函数 6 epsilonT=1e-12;%收敛判断标准：相对误差 7 x0=pi/2;%给初值 8 x1=pi/2+0.1; 9 n=0;%迭代次数 10 while 1 11 % x2=x0-f(x0)*(x1-x0)/(f(x1)-f(x0)); 12 x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0));%等同于上句，都行 13 epsilon=abs((x1-x0)/x1); 14 x0=x1; 15 x1=x2; 16 n=n+1; 17 disp([\'第\' num2str(n) \'次迭代，x=\', num2str(x1,15)]); 18 if epsilon<epsilonT 19 break;%跳出while循环 20 end 21 end 22 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
非线性多元方程组用牛顿法求解：
1 % https://zhuanlan.zhihu.com/p/146371408 2 % https://zhuanlan.zhihu.com/p/63103354 知乎代码 3 clc; 4 clear all; 5 close all; 6 x0=[1 2]; 7 eps=1e-12; 8 [allx,ally,r,n]=mulNewton(fun,x0,eps); 9 10 disp([\'迭代\' num2str(n) \'次，\' \'x1=\' num2str(eval(r(1)),100) \',x2=\' num2str(eval(r(2)),100) ])% 0.19808577588668504464406945200284 11 % disp(r) 12 disp(\'newton end\') 13 %% 子函数区域 14 function [allx,ally,r,n]=mulNewton(F,x0,eps) 15 if nargin==2 16 eps=1.0e-4; 17 end 18 x0 = transpose(x0); 19 Fx = subs(F,transpose(symvar(F)),x0);%将初始点代入方程组 20 var = transpose(symvar(F)); 21 dF = jacobian(F,var);%求雅克比矩阵 22 dFx = subs(dF,transpose(symvar(F)),x0);%将当前解带入雅克比矩阵 23 r=x0-inv(dFx)*Fx\';%迭代 24 n=1; 25 tol=1; 26 N=100; 27 symx=length(x0); 28 ally=zeros(symx,N); 29 allx=zeros(symx,N); 30 31 while tol>eps 32 x0=r; 33 Fx = subs(F,transpose(symvar(F)),x0); 34 dFx = subs(dF,transpose(symvar(F)),x0); 35 r=vpa(x0-inv(dFx)*Fx\'); 36 tol=norm(r-x0) 37 if(n>N) 38 disp(\'迭代步数太多，可能不收敛！\'); 39 break; 40 end 41 allx(:,n)=x0; 42 ally(:,n)=Fx; 43 n=n+1; 44 end 45 end 46 47 48 function f = fun(x) 49 k=2; 50 for i=1:k 51 x(i)=sym ([\'x\',num2str(i)]); 52 end 53 54 f(1)=0.5*sin(x(1))+0.1*cos(x(1)*x(2))-x(1); 55 f(2)=0.5*cos(x(1))-0.1*cos(x(2))-x(2); 56 end
1stopt求解非线性方程组：
1 NewCodeBlock"求解非线性方程组"; 2 //https://zhuanlan.zhihu.com/p/63103354 3 //https://zhuanlan.zhihu.com/p/146371408 4 AlgorithmOption =[1,0,1.00E-30,1000,20,30,50,20,0,0.15,1]; 5 Parameter x(2); 6 Function 0.5*sin(x1)+0.1*cos(x1*x2)-x1=0; 7 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
高斯列主元消去法：
1 % 高斯列主元消去法 https://zhuanlan.zhihu.com/p/128577559 2 % 先选择每列元素绝对值最大值放通过行变换放在主对角线上，之后将矩阵A化为上三角矩阵，然后回代求解线性方程组Ax=b。 3 clc; 4 clear all; 5 close all; 6 % A=[2 8 2; 7 % 1 6 -1; 8 % 2 1 2];%https://wenku.baidu.com/view/99bf58cf050876323112120c.html 9 % b=[14 13 5]\'; 10 A=[0.001 2 3; 11 -1 3.712 4.623; 12 -2 1.072 5.643];%https://wenku.baidu.com/view/99bf58cf050876323112120c.html 13 b=[1 14 2 15 3 ]; 16 x = gaussPivotScale(A, b) 17 disp(\'end\'); 18 19 function x = gaussPivotScale(A, b) 20 21 %GAUSSPIVOTSCALE It uses Gauss elimination with partial pivoting 22 % and row scaling to solve the linear system Ax = b. 23 % A : It is the n-by-n coefficient matrix. 24 % b : It is the n-by-1 result vector. 25 % x : It is the n-by-1 solution vector. 26 27 tic 28 format long 29 30 % judge wether there is a solution or not from the linear system 31 r_A = rank(A); 32 r_Ab = rank([A, b]); 33 34 if r_A ~= r_Ab 35 disp(\'This linear system is no solution\'); 36 end 37 38 % Elimination 39 40 [row, ~] = size(A); 41 C = [A, b]; 42 43 for k = 1 : row-1 44 45 % style of table 46 fprintf(\'This is%2dth transformation \n [A | b] = \n\', k); 47 48 M = max(abs(C(k:end, k:end-1)), [], 2); % find maximum in kth row,A所有行每行元素绝对值最大值 49 a = abs(C(k:end, k)); % each value is absoluted value in kth column，所有行第k列绝对值 50 [~, id] = max(a ./ M); % find row with maximum 51 52 if id > k 53 C([k, id], :) = C([id, k], :); % pivot rows 54 end 55 56 mult = C(k+1:row, k) / C(k, k); % construct multipliers 57 58 % creat mesh, and improve speed 59 [C_k, mult] = meshgrid(C(k, :), mult); 60 C(k+1:row, :) = C(k+1:row, :) - C_k .* mult; 61 62 disp(C) 63 end 64 65 % Back Substitution 66 [row, ~] = size(C); 67 68 for k = row : -1 : 1 69 C(k, :) = C(k, :) ./ C(k, k); 70 C(1:k-1, row+1) = C(1:k-1, row+1) - C(1:k-1, k) * C(k, row+1); 71 end 72 73 x = C(:, end); 74 75 toc 76 77 end
