第二章 非线性方程（组）的数值解法

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 第二章 非线性方程（组）的数值解法 2.1 方程（组）的根及其MATLAB命令

2.1.2 求解方程（组）的solve命令

求方程f(x)=q(x)的根可以用MATLAB命令:

>> x=solve('方程f(x)=q(x)','待求符号变量x')

求方程组fi(x1,…,xn)=qi(x1,…,xn) (i=1,2,…,n)的根可以用MATLAB命令:

>>E1=sym('方程f1(x1,…,xn)=q1(x1,…,xn)'); …………………………………………………….

En=sym('方程fn(x1,…,xn)=qn(x1,…,xn)'); [x1,x2,…,xn]=solve(E1,E2,…,En, x1,…,xn)

2.1.4 求解方程（组）的fsolve命令

fsolve的调用格式： X=fsolve(F,X0)

2.2 搜索根的方法及其MATLAB程序

2.2.1 作图法及其MATLAB程序

作函数y?f(x)在区间 [a,b] 的图形的MATLAB程序一

x=a:h:b; % h是步长 y=f(x); plot(x,y)

grid, gtext('y=f(x)')

说明：⑴ 此程序在MATLAB的工作区输入，运行后即可出现函数y?f(x)的图形.此图形与x轴交点的横坐标即为所要求的根的近似值.

⑵ 区间[a,b] 的两个端点的距离 b-a 和步长h的绝对值越小，图形越精确. 作函数y?f(x)在区间 [a,b]上的图形的MATLAB程序二 将y?f(x)化为h(x)?g(x)，其中h(x)和g(x)是两个相等的简单函数.

x=a:h:b; y1=h(x); y2=g(x); plot(x, y1, x, y2)

grid,gtext(' y1=h(x),y2=g(x)') 说明：此程序在MATLAB的工作区输入，运行后即可出现函数y1?h(x)和y2?g(x)的图形.两图形交点的横坐标即为所要求的根的近似值.

2.2.2 逐步搜索法及其MATLAB程序

逐步搜索法的MATLAB主程序

function [k,r]=zhubuss(a,b,h,tol)

% 输入的量--- a和b是闭区间[a,b]的左、右端点； %---h是步长；

%---tol是预先给定的精度.

% 运行后输出的量---k是搜索点的个数；

7．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 % --- r是方程 在[a,b]上的实根的近似值，其精度是tol； X=a:h:b;Y=funs(X);n=(b-a)/h+1;m=0; X(n+1)=X(n);Y(n+1)=Y(n); for k=2:n

X(k)=a+k*h;Y(k)=funs(X(k)); %程序中调用的funs.m为函数 sk=Y(k)*Y(k-1); if sk<=0,

m=m+1;r(m)=X(k); end

xielv=(Y(k+1)-Y(k))*(Y(k)-Y(k-1)); if (abs(Y(k))

例2.2.4 用逐步搜索法的MATLAB程序分别求方程2x?2x32?3x?3?0和

sin(cos2x3)?0在区间[?2,2]上的根的近似值，要求精度是0.000 1.

解 将逐步搜索法的MATLAB程序保存名为zhubuss.m的M文件. 建立M文件funs.m

function y=funs(x)

y=2.*x.^3+2.*x.^2-3.*x-3 在MATLAB工作窗口输入如下程序

>> [k,r]=zhubuss(-2,2,0.001,0.0001)

运行后输出的结果

k =4001

r = -1.2240 -1.0000 -1.0000 -0.9990 1.2250

即搜索点的个数为k =4 001，其中有5个是方程2x?2x?3x?3?0的近似根，即r = -1.224 0，-1.000 0，-1.000 0，-0.999 0，1.225 0，其精度为0.000 1.

在程序中将y=2.*x.^3+2.*x.^2-3.*x-3用y=sin(cos(2.*x.^3)) 代替，可得到方程

32sin(cos2x3)?0在区间[?2,2]上的根的近似值如下

r = -1.9190 -1.7640 -1.5770 -1.3300 -0.9220 0.9230

1.3310 1.5780 1.7650 1.9200

2.3 二分法及其MATLAB程序

2.3.2 二分法的MATLAB程序 二分法的MATLAB主程序

function [k,x,wuca,yx]=erfen(a,b,abtol) a(1)=a; b(1)=b;

ya=fun(a(1)); yb=fun(b(1)); %程序中调用的fun.m 为函数 if ya* yb>0,

disp('注意：ya*yb>0,请重新调整区间端点a和b.'), return end

max1=-1+ceil((log(b-a)- log(abtol))/ log(2)); % ceil是向?? 方向取整

for k=1: max1+1

a;ya=fun(a); b;yb=fun(b); x=(a+b)/2; yx=fun(x); wuca=abs(b-a)/2; k=k-1; [k,a,b,x,wuca,ya,yb,yx] if yx==0

a=x; b=x; elseif yb*yx>0

b=x;yb=yx; else

8．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 a=x; ya=yx; end

if b-a< abtol , return, end end

k=max1; x; wuca; yx=fun(x);

例2.3.3 确定方程x3-x+4=0的实根的分布情况，并用二分法求在开区间 (-2,-1)内的实根的近似值，要求精度为0.001.

解 （1）先用两种方法确定方程x3-x+4=0 的实根的分布情况。 方法1 作图法.

在MATLAB工作窗口输入如下程序

>>x=-4:0.1:4;

y=x.^3-x +4; plot(x,y) grid,gtext('y=x^3-x+4') 画出函数f(x)=x3-x+4的图像.从图像可以看出，此曲线有两个驻点?3都在x轴的上方，3在(-2,-1)内曲线与x轴只有一个交点，则该方程有唯一一个实根，且在 (-2,-1)内.

方法2 试算法.

在MATLAB工作窗口输入程序

>>x=-4: 1:4,y=x.^3-x +4

运行后输出结果

x = -4 -3 -2 -1 0 1 2 3 4 y = -56 -20 -2 4 4 4 10 28 64

由于连续函数f(x)满足f(?2)?f(?1)?0，所以此方程在(-2,-1)内有一个实根.

（2） 用二分法的主程序计算. 在MATLAB工作窗口输入程序

>> [k,x,wuca,yx]=erfen (-2,-1,0.001)

运行后屏幕显示用二分法计算过程被列入表 2-3，其余结果为

k = 9，x=-1.7959，

wuca = 9.7656e-004，yx = 0.0037

表 2-3 次数k 左端点ak 右端点bk 0 1 2 3 4 5 6 7 8 9 -2.000 0 -2.000 0 -2.000 0 -1.875 0 -1.812 5 -1.812 5 -1.796 9 -1.796 9 -1.796 9 -1.796 9 -1.000 0 -1.500 0 -1.750 0 -1.750 0 -1.750 0 -1.781 3 -1.781 3 -1.789 1 -1.793 0 -1.794 9 中点xk -1.500 0 -1.750 0 -1.875 0 -1.812 5 -1.781 3 -1.796 9 -1.789 1 -1.793 0 -1.794 9 -1.795 9 bk?ak 函数值2f(ak) 0.500 0 0.250 0 0.125 0 0.062 5 0.031 3 0.015 6 0.007 8 0.003 9 0.002 0 0.001 0 函数值f(bk) -2.000 0 4.000 0 -2.000 0 -2.000 0 -0.716 8 -0.141 8 -0.141 8 -0.004 8 -0.004 8 -0.004 8 -0.004 8 2.125 0 0.390 6 0.390 6 0.390 6 0.129 6 0.129 6 0.062 7 0.029 0 0.012 1 函数值f(xk) 2.125 0 0.390 6 -0.716 8 -0.141 8 0.129 6 -0.004 8 0.062 7 0.029 0 0.012 1 0.003 7 2.4 迭代法及其MATLAB程序

2.4.2 迭代法的MATLAB程序1

迭代法的MATLAB主程序1

function [k,piancha,xdpiancha,xk]=diedai1(x0,k) % 输入的量--x0是初始值,k是迭代次数

9．

第三篇 第二章 非线性方程（组）的数值解法的MATLAB程序 x(1)=x0; for i=1:k

x(i+1)=fun1(x(i));%程序中调用的fun1.m为函数y=φ(x) piancha= abs(x(i+1)-x(i));

xdpiancha=piancha/( abs(x(i+1))+eps);

i=i+1;xk=x(i);[(i-1) piancha xdpiancha xk] end

if (piancha >1)&(xdpiancha>0.5)&(k>3)

disp('请用户注意：此迭代序列发散，请重新输入新的迭代公式') return; end

if (piancha < 0.001)&(xdpiancha< 0.0000005)&(k>3) disp('祝贺您！此迭代序列收敛，且收敛速度较快') return; end

p=[(i-1) piancha xdpiancha xk]';

例2.4.1 求方程f(x)?x2?2x?10的一个正根. 解 在MATLAB工作窗口输入程序

>> [k,piancha,xdpiancha,xk]= diedai1(2,5)

2运行后输出用迭代公式xk?1?(10?xk)/2的结果

[k,piancha,xdpiancha,xk]=

1.00000000000000 1.00000000000000 0.33333333333333 3.00000000000000

2.00000000000000 2.50000000000000 5.00000000000000 0.50000000000000

3.00000000000000 4.37500000000000 0.89743589743590 4.87500000000000

4.00000000000000 11.75781250000000 1.70828603859251 -6.88281250000000

5.00000000000000 11.80374145507813 0.63167031671297 -18.68655395507813

请用户注意：此迭代序列发散，请重新输入新的迭代公式

k=5，piancha = 11.80374145507813，xdpiancha = 0.63167031671297， xk = -18.68655395507813

*由以上运行后输出的迭代序列与根x=2.316 624 790 355 40相差越来越大，即迭代序列{xk}发散，此迭代法就失败.这时偏差piancha逐渐增大且偏差的相对误差xdpiancha的值大于0.5.

用迭代公式xk?1?10/(xk?2)运行后输出的结果 [k,piancha,xdpiancha,xk]=

1.00000000000000 0.50000000000000 0.20000000000000 2.50000000000000

2.00000000000000 0.27777777777778 0.12500000000000 2.22222222222222

3.00000000000000 0.14619883040936 0.06172839506173 2.36842105263158

4.00000000000000 0.07926442612555 0.03462603878116 2.28915662650602

5.00000000000000 0.04230404765128 0.01814486863115 2.33146067415730

k =5,piancha =0.04230404765128,xdpiancha = 0.01814486863115, xk = 2.33146067415730

可见，偏差piancha和偏差的相对误差xdpiancha的值逐渐变小，且第5次的迭代值xk =2.331 460 674 157 30与根x=2.316 624 790 355 40接近，则迭代序列{xk}收敛，但收敛速度较慢，此迭代法较为成功.

用迭代公式xk?1?xk?(xk?2xk?10)/(2xk?2)运行后输出的结果

10．

*2