第七讲 MATLAB中求方程的近似根()
教学目的:学习matlab中求根命令,了解代数方程求根求解的四种方法,即图解法、准解析法、数值方法以及迭代方法,掌握对分法、迭代法、牛顿切法线求方程近似根的基本过程;掌握求代数方程(组)的解的求解命令.
教学重点:求方程近似解的几种迭代方法代数方程(组)的解的求解命令的使用方法.利用所学的编程知识,结合具体的实例,编制程序进行近似求根.掌握相关的代数方程(组)的求解命令及使用技巧.
教学难点:方程的近似求解和非线性方程(组)的求解.
一、问题背景和实验目的
求代数方程的根是最常见的数学问题之一(这里称为代数方程,主要是想和后面的微分方程区别开.为简明起见,在本实验的以下叙述中,把代数方程简称为方程),当是一次多项式时,称为线性方程,否则称之为非线性方程.
是非线性方程时,由于的多样性,尚无一般的解析解法可使用,但如果对任意的精度要求,能求出方程的近似根,则可以认为求根的计算问题已经解决,至少能满足实际要求.同时对于多未知量非线性方程()而言,简单的迭代法也是可以做出来的,但在这里我们介绍相关的命令来求解,不用迭代方法求解.
通过本实验,达到下面目的:
1. 了解对分法、迭代法、牛顿切线法求方程近似根的基本过程;
2. 求代数方程(组)的解.
首先,我们先介绍几种近似求根有关的方法
1. 对分法
对分法思想:将区域不断对分,判断根在某个分段内,再对该段对分,依此类推,直到满足精度为止.对分法适用于求有根区间内的单实根或奇重实根.
上连续,,即.则根据连续函数的介值定理,在内至少存在一点,使
下面的方法可以求出该根:
(1) ,计算
(2) ,则的根,停止计算,输出结果
,则令,若,则令.……,有以及相应的
(3) (为预先给定的精度要求),退出计算,输出结果
反之,返回(1),重复(1)(2)(3)
以上方法可得到每次缩小一半的区间序列,在中含有方程的根.
当区间长很小时,取其中点为根的近似值,显然有
以上公式可用于估计对分次数
分析以上过程不难知道,对分法的收敛速度与公比为的等比级数相同.由于,可知大约对分10次,近似根的精度可提高三位小数.对分法的收敛速度较慢,它常用来试探实根的分布区间,或求根的近似值.
2. 迭代法
a) 松弛法:由方程构造一个等价方程
则迭代方程是:
,其中
松弛法的加速效果是明显的 (见附录4),甚至不收敛的迭代函数经加速后也能获得收敛.
b) Altken方法:
松弛法要先计算,在使用中有时不方便,为此发展出以下的 Altken 公式:
这就是Altken 公式,它的加速效果也是十分明显的,它同样可使不收敛的迭代格式获得收敛(见附录5)
3. 牛顿(Newton)法(牛顿切线法)
是非线性方程
其迭代公式为:
 
即为牛顿法公式.
牛顿法的缺点是:(1)对重根收敛很慢;(2)对初值要求较严,要求相当接近真值
因此,常用其他方法确定初值,再用牛顿法提高精度.
以下是本实验中的几个具体的实验
具体实验1:对分法
先作图观察方程:的实根的分布区间,再利用对分法在这些区间上分别求出根的近似值.程序如下:
function [y,p]=erfen()
clc, x=[];a=[];b=[]; a(1)=1;b(1)=2; i=1;x(i)=(a(i)+b(i))/2; e=abs(f(x(i)));
ezplot('x^3-3*x+1',[a(1),b(1)]);hold on,   plot([a(i),b(i)],[0,0])
while e>10^(-5)   
    plot([a(i),a(i)],[0,100],[x(i) x(i)],[0 100],[b(i) b(i)],[0 100]),pause(0.5)
    if f(a(i))*f(x(i))<0
        a(i+1)=a(i);b(i+1)=x(i);x(i+1)=(a(i+1)+b(i+1))/2;
    else
        a(i+1)=x(i);b(i+1)=b(i);x(i+1)=(a(i+1)+b(i+1))/2;
    end
    e=abs(f(x(i)));i=i+1;
end
y=x(i);p=[a;x;b]'
function u=f(x)
u=x^3-3*x+1;
end
end
图形如下:
结果为:1.5321
具体实验2:普通迭代法
采用迭代过程:求方程 0.5 附近的根,精确到第 4 位小数.
构造等价方程:
用迭代公式:
具体实验3:迭代法的加速1——松弛迭代法
迭代公式为 
clc;x=[];w=[]; x(1)=1;w(1)=1/(1-x(1));
for i=1:10
    w(i)=1/(1- x(i));    x(i+1)=(1-w(i))*x(i)+ w(i)*(x(i)^3+1)/3;
end
x
另外有程序可以参考,详见参见附录4
具体实验4:迭代法的加速2——Altken迭代法
迭代公式为:
%(符号计算)
syms x fx gx;
gx=(x^3+1)/3;fx=x^3-3*x+1; disp('k    x    x1    x2')
x=0.5;k=0; ffx=subs(fx, 'x', x);
while abs(ffx)>0.0001;
u=subs(gx, 'x', x);v=subs(gx, 'x', u);
  disp([num2str(k), '  ', num2str(x), '  ', num2str(u), '  ', num2str(v)])
x=v-(v-u)^2/(v-2*u+x);k=k+1;ffx=subs(fx, 'x', x);
end
disp([num2str(k), '  ', num2str(x), '  ', num2str(u), '  ', num2str(v)])
%(数值计算)
function [y,p]=althken() 
% 求方根的迭代程序
clc,format long e, x(1)=6; i=1;p=[];ezplot('x^3-3*x+1',[x(1)-9,x(1)+1]);hold on
  plot([x(1)-20,x(1)+2],[0,0])
while abs(f(x(i)))>=10^(-5) 
      plot(x(i),0,'*')
      t1=phi(x(i));t2=phi(t1);      x(i+1)=t2-(t2-t1)^2/(t2-2*t1+x(i)+eps);
      p=[p;[i, x(i),t1,t2]];  i=i+1;      pause(0.1)
end 
p,y=x(i), i, format
function u=phi(x)
    u=(x^3+1)/3;
end
  function u=f(x)
    u=x^3+1-3*x;
  end
end
具体实验5:牛顿法
用牛顿法计算方程-22之间的三个根.
提示:迭代公式:
function [y,p]=newton() 
% 求方根的迭代程序
clc,format long e, x(1)=6; i=1; p=[]; ezplot('x^3-3*x+1',[x(1)-9,x(1)+1]);hold on
  plot([x(1)-20,x(1)+2],[0,0])
while abs(f(x(i)))>=10^(-5) 
      plot(x(i),0,'*'),  x(i+1)=x(i)-f(x(i))/(df(x(i))+eps);
      p=[p;[i, x(i)]];  i=i+1;      pause(0.1)
end 
format short, p,y=x(i), i, 
function u=df(x)
  u=3*x^2-3;
end
function u=f(x)
  u=x^3+1-3*x;
  end
end
结果:
结果为:  1.5321
   
※进一步思考:用迭代法求3的平方根. 迭代公式为. 编写M函数文件My_sqrt.m, 3正的平方根. 要求误差小于.仅要求写出源程序.试使用以上介绍的迭代法来相互比较
参考程序:
function y=my_sqrt(a)                                                   
% 求方根的迭代程序
if nargin~=1|~isa(a,'double') ,  error('输入数字为一个正数!'),end
if a<0,  error('输入数字为正数!'), end
if a>0                                                     
  format long e,  x(1)=0;  x(2)=1; i=1;                                   
  while abs(x(i+1)-x(i))>=10^(-5)
matlab学好了有什么用i=i+1;x(i+1)=1/2*(x(i)+a/(x(i)+eps));end 
      y=x(i+1);i,format                                                       
end
现在我们简单介绍图解法如何来求解一元方程和二元方程的根
例:exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)=0.5
>>ezplot('exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5',[0 5])
>>hold on, line([0,5],[0,0])

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。