数值计算大作业:非线性方程求根(二分法、牛顿法、弦截法在Matlab实现)

发布于:2022-11-28 ⋅ 阅读:(317) ⋅ 点赞:(0)

作为研究生的入门课,数值计算的大作业算是所有研究生开学的重要编程作业。

    我把二分法、牛顿法、弦截法求解非线性方程求根的数值计算作业在MATLAB中编程实现。具体的程序详细标注后放在文章附录了,算法数学原理也一并放在文中了。需要的同学自取。下文为作业详解

题目:给定方程:sin10x-xe-x=x,完成以下工作

1.绘出函数f(x)=sin10x-xe-x-x的图形,确定有根区间。

函数图像如下图,可以确定函数有三个根,区间分别为[-0.3,-0.2],[-0.1,0.1],[0.2,0.3].

2.使用二分法求出方程的根,要求精度ε10-8

二分法数学原理:

对于区间[a,b]上连续不断且f(a)*f(b)<0的函数y=f(x),通过不断地把函数f(x)的零点所在的区间一分为二,使区间的两个端点逐步逼近零点,进而得到零点近似值的方法叫二分法。

给定精确度ξ,用二分法求函数f(x)零点近似值的步骤如下:

1 确定区间[a,b],验证f(a)·f(b)<0,给定精确度ξ.

2 求区间(a,b)的中点c.

3 计算f(c).

(1) 若f(c)=0,则c就是函数的零点;

(2) 若f(a)·f(c)<0,则令b=c;

(3) 若f(c)·f(b)<0,则令a=c.使用 Newton 法求出方程的根,要求精度

(4) 判断是否达到精确度ξ:即若|a-b|<ξ,则得到零点近似值a(或b),否则重复2-4.

具体运行程序在附录中,此处只展示结果

计算结果如下表所示

有根区间

二分次数

f(x)=0的根x*

f(x*)

(-0.3,-0.2)

28

-0.2526

3.4673*10-4

(-0.1,0.1)

28

0

0

(0.2,0.3)

28

0.2654

-4.3458*10-4

 3.使用Newton法求出方程的根,要求精度ε10-8

牛顿法数学原理:

     对求解方程f(x)=0,即f(x0)+(x-x0)*f'(x0)=0,求解x:x1=x0-f(x0)/f'(x0),因为这是利用泰勒公式的一阶展开,f(x)=f(x0)+(x-x0)f'(x0)处并不是完全相等,而是近似相等,这里求得的x1并不能让f(x)=0,只能说f(x1)的值比f(x0)更接近f(x)=0,所以,迭代求解的想法可以实现。进而推出x(n+1)=x(n)-f(x(n))/f'(x(n)),通过迭代,这个式子必然在f(x*)=0的时候收敛

计算结果如下表所示

迭代初值

迭代次数

f(x)=0的根x*

f(x*)

-0.3

4

-0.2526

3.4673*10-4

-0.5

9

1.3208*10-8

-1.0566*10-8

0.2

14

0.2654

-4.3458*10-4

4.使用弦截法求出方程的根,要求精度ε10-8

  弦截法数学原理:

    弦截法是求非线性方程近似根的一种线性近似方法。它是以与曲线弧AB对应的弦AB与x轴的交点横坐标作为曲线弧AB与x轴的交点横坐标的近似值μ来求出方程的近似解。该方法一般通过计算机编程来实现。弦截法的原理是以直代曲即用弦(直线)代替曲线求方程的近似解,也就是利用对应的弦与轴的交点横坐标来作为曲线弧与轴的交点横坐标的近似值。

   Newton法要计算函数的导数,当导数不方便计算时,可以利用导数的定义,由相近点处函数值的差来近似,这就得到弦截法公式:

 

弦截法原理图

 

计算结果如下表所示

迭代初值

迭代次数

f(x)=0的根x*

f(x*)

(-0.3,-0.2)

7

-0.2526

3.4673*10-4

(-0.1,0.1)

10

1.3208*10-8

-1.0566*10-8

(0.2,0.3)

28

0.2654

-4.3458*10-4

结论:

     三种非线性方程求根方法只要满足各自的收敛条件,最后求得的根都很接近。但是牛顿法迭代的次数和是否收敛过度依赖迭代初值的选择,而弦截法即避免了初值的选取,也避开了求原函数的导数,大大减少了计算量。

5.选择方程的一个根,求出以上3种方法的数值收敛阶。(选做题)

PS:我字太丑了,但是结果是对的 

附录:Matlab程序(分为主体程序,二分法函数程序,牛顿法函数程序,弦截法函数程序)

主体程序

% 非线性方程求根
%Solution of Non-linear Equation
% 绘制函数图像
f=@(x) sin(10*x)-x.*exp(-x)-x;
x=-0.5:0.01:0.5;
y=f(x);
plot(x,y)
grid on;
% 多次画图,确定此函数有三个根,因函数非单调,划分为3个有根区间[-0.3,-0.2],[-0.1,0.1],[0.2,0.3]
% 用二分法求这三个区间的方程解
bis1=my_bisection(-0.3,-0.2);
bis2=my_bisection(-0.1,0.1);
bis3=my_bisection(0.2,0.3);
% 用牛顿法求这三个区间的方程解
[root_new1,kn1] = newton(-0.2);
[root_new2,kn2] = newton(-0.25);
[root_new3,kn3] = newton(0.3);
% 用弦截法求这三个区间的方程解
[root_sec1,ks1] = secant_method(-0.3,-0.2);
[root_sec2,ks2] = secant_method(-0.1,0.1);
[root_sec3,ks3] = secant_method(0.2,0.3);

二分法函数程序

function root = my_bisection(a,b)
%root为二分法输出的根,输入的a,b为判断出的有根区间
% 误差精度要低于10^(-8);
accuracy=1.0e-8;
% 本题求根函数
f=@(x) sin(10*x)-x.*exp(-x)-x;
% 计算左端点函数值
fa=f(a);
% 计算右端点函数值
fb=f(b);
% 计算中点函数值
fmid=f((a+b)/2);

% 判断函数情况,改变函数计算区间
% 使用递归方法
if(fa*fmid>0) 
    a_bis=(a+b)/2;
    root=my_bisection(a_bis,b);
else
    if(fa*fmid==0)
      root=(a+b)/2;
    else
     if(abs(b-a)<=accuracy)
        root=(a+b)/2;
     else
       b_bis=(a+b)/2;  
       root=my_bisection(a,b_bis);
     end
    end
end

牛顿法函数程序

function [root,k] = newton(x0)
%root为牛顿迭代法输出的根,k为迭代次数,输入的x0为迭代初值
% 误差精度要低于10^(-8);
accuracy=1.0e-8;
% 原函数和相应的一阶导数
syms x
f=sin(10*x)-x.*exp(-x)-x;
df=diff(sym(fx));
% ddf=diff(sym(df));
% f='sin(10*x)-x.*exp(-x)-x';
% df='10*cos(10*x)-exp(-x)+x.*exp(-x)-1';
% 初值准备,其中误差项error_term为1是为了之后进去while循环的初始条件
x1=x0;
k=0;
error_term=0.1;
        while((error_term>accuracy)&&(k<=300))
            k=k+1;
            fx=subs(fx,x1);
            df=subs(df,x1);
            root=x1-fx/df;
            error_term=abs(root-x1);
            x1=root;
        end
        
end

弦截法函数程序

function [root,k] = secant_method(a,b)
%root为弦截法输出的根,k为迭代次数,输入的a,b为区间的两个节点
% 误差精度要低于10^(-8);
accuracy=1.0e-8;
% 本题求根函数
f=@(x) sin(10*x)-x.*exp(-x)-x;
% 计算左端点函数值
fa=f(a);
% 计算右端点函数值
fb=f(b);
%每一步弦截的跟
root=a-(b-a)*fa/(fb-fa);
% 初值准备,其中误差项error_term为1是为了之后进去while循环的初始条件
k=0;
error_term=0.1;
        while((error_term>accuracy)&&(k<=300))
            k=k+1;
           x1=root;
           fx=f(root);
%确定衰减方向,保证是往f减少出缩进
            s=fx*fa;
           if(s>0)
               root=b-(x1-b)*fb/(fx-fb);
           else
               root=a-(x1-a)*fa/(fx-fa);
           end
            error_term=abs(root-x1);
        end
end

本文含有隐藏内容,请 开通VIP 后查看

网站公告

今日签到

点亮在社区的每一天
去签到