解非线性方程

实验类型：●验证性实验  ○综合性实验  ○设计性实验

实验目的：进一步熟练掌握解非线性方程的二分法算法、牛顿迭代法，提高编程能力和解算非线性方程问题的实践技能。

实验内容：用二分法算法(取[a,b]=[1,2])、牛顿迭代法解算非线性方程

的根（迭代计算时要求取精度ε=10-5）

实验原理   二分法算法

           牛顿迭代法

实验步骤

1）要求上机实验前严格按照二分法的算法步骤对于非线性方程

手工迭代计算三步，并把计算过程步骤详细写在A4纸上；对手工计算做进一步的总结梳理然后画出二分法程序设计流程图；应用Matlab语言在纸上编写出二分法的程序代码。

2）要求上机实验前严格按照牛顿迭代算法步骤对于非线性方程

3） 编辑录入程序

4）输入数据（截图），调试程序并记录调试过程中出现的问题（截图）及修改程序的过程（截图），把你的截图及时粘贴在Word文档中，并对截图配上必要的文字说明，以便证明和描述你的实验过程真实可靠。

5）经反复调试后，运行程序并验证程序运行的前三步中间结果和手工计算的前三步中间结果是否一致，要求通过列表对比、分析和论证程序运行的前三步中间结果和手工计算的前三步中间结果是否一致，若不一致，必须找到不一致的原因，并预以纠正，最后明确给出对比分析结论。 

实验思考

1. 试验改变牛顿迭代算法的迭代初值进行计算，分别对比计算结果并作出分析总结
2. 试用实验证实结论：在计算机（如Matlab上）上实际的的浮点算术体系中二分法运行结果的准确度不可能随迭代过程一直提高。即要用实验证实：在精度ε小于等于Matlab的固定变量eps的情况下，二分法最终区间[a,b]的两端点是两个相邻的浮点数。

数学建模与应用：请建立数学模型并应用你的程序解决一个物理计算问题--木质球体浸入水中深度问题。假若球体是由密度为ρ=0.638且半径为r=10cm的长叶松木质球体构成的，问球体浸入水中深度是多少？又假若球体是由密度为ρ=0.638且三个长半轴分别为a=8cm,b=10cm,c=12cm的长叶松木质椭球体构成的，问椭球体浸入水中深度是多少？【提示：1）通过受力分析建立力学平衡方程--非线性方程，用二分法和牛顿迭代法解算非线性方程，阐释计算结果的物理意义；2）注意椭球体的在水中可能有多种平衡方式，在水中的平衡方式不同，其数学模型可能不同，相应的解答也会有差异，请针对各种情况建模并解答。迭代计算时要求取精度ε=10-5。】

实验报告书写要求：根据实验情况和结果撰写并递交实验报告。实验报告应当有算法原理简介，算法流程图，程序代码，运行调试记录，计算结果列表对比分析，对比分析的结论，从物理角度对计算结果做出阐释；对于数学建模的应用问题要有分析建模的详细过程、求解计算过程、实际问题解答或阐释的内容。在A4纸上对每个问题分别用二分法和牛顿迭代法手工迭代计算前三步，手工计算要求所有数据保留5位小数，然后将手工迭代计算每步的计算结果与计算机程序计算的前三步每步结果进行比较，若数据整数及小数点后4位小数都相同则认为是相同的，若不同，要分析不同的原因并予以校正，经分析、校正后给出对比验证的结论。 将手工计算过程及其验证结论拍照成JPG图片插入在实验报告的word 电子版中。

实验报告打印和装订顺序要求

实验报告打印要求：在A4纸上将实验报告的Word文档双面打印；打印的字迹要清晰，报告文字或截图上的文字打印后不得小于宋体6号，要确保字迹清晰可见。

实验报告装订顺序

• 实验任务（本文档版式不得做任何改动）；
• 实验报告正文（Word文档）；
• 可以将A4纸上详细手工计算过程拍照成JPG图片插入在实验报告的word 电子版中随Word文档双面打印。
• 左边缘装订

实验总结(学会了......; 掌握了......; 训练了......; 发现了......; 今后学习中......有待提高。）

附件：程序代码

电子报告word文件命名规则：专业班级-完整学号-实验X-姓名.doc, 如信息123班学号为201212030315的郭海涛同学实验2报告word文件命名则应是：信息123-201212030315-实验2-郭海涛.doc, 其中 .doc是Word文件扩展名。

特别提醒：学委负责检查电子文件命名规范，命名不规的将不接收。每人必需提交电子报告和纸质报告各一份(两者内容一致)。本实验任务书版式不得做任何改动。电子报告以班为单位提交压缩包.

此页之前的内容为实验任务内容，实验任务内容及其版式不得做任何更改；实验报告从下页开始编写排版。

手工计算迭代与Matlab迭代对比

手工计算结果保留五位小数，计算机结果保留四位小数，经过比对，结果完全一致。

• 用牛顿迭代法解算非线性方程的根

程序录入及运行界面截图：

手工计算迭代与Matlab迭代对比

手工计算结果保留五位小数，计算机结果保留四位小数，经过比对，结果完全一致。

• 实验思考

1.试验改变牛顿迭代算法的迭代初值进行计算，分别对比计算结果并作出分析总结

我改变了初始迭代值分别为1，1.5，2迭代结果如图所示，会有不同，但最后还是会在1.21341这个解这里会合。迭代次数也会有所改变，当初始值变为2时会比其他两个初值多迭代一次。（还有可能改变迭代收敛性）

2.

a = 1;b = 2;tol = eps;  maxK= 100;

[root, k, a, b] = bi(a, b, tol, maxK);

fprintf('近似根为: %f\n', root);fprintf('对应的函数值为: %f\n', f(root));

fprintf('最终区间的左端点 a = %0.20f\n', a);fprintf('最终区间的右端点 b = %0.20f\n', b);

function [root, k, a, b] = bi(a, b, tol, maxK)

    k = 0;

    while k < maxK

        c = (a + b) / 2;fc = f(c);

        if abs(fc) < tol

            root = c;return;

        elseif f(a) * fc < 0

            b = c;

        else

            a = c;

        end

        k = k + 1;

    end

    root = (a + b) / 2;

end

function y = f(x)

    y = x^3 + x - 3;

end

由实验结果可以印证：在计算机（如Matlab上）上实际的的浮点算术体系中二分法运行结果的准确度不可能随迭代过程一直提高。即在精度ε小于等于Matlab的固定变量eps的情况下，二分法最终区间[a,b]的两端点是两个相邻的浮点数。

• 数学建模

模型建立：

F浮=G球；

即：

球体：(单位：m)

则有：f(x)=-x^3+0.3x^2-0.002552,x∈[0,0.2]

对应二分法和牛顿迭代法结果如下：

椭球体：

当c长半轴与水平面平行时：（单位：cm）

由公式代入数值有：f(x)= 25x^3-120x^2-6717,x∈[0,35]

当b长半轴与水平面平行时：（单位：cm）

f(x)=8*x^3-120*x^2-10183,x∈[0,35]

当a长半轴与水平面平行时：（单位：cm）

f(x)=5*x^3-36*x^2-4399,x∈[0,35]

• 实验总结

学会了二分法迭代和牛顿迭代求解非线性方程组，以及这两种解法的的具体迭代过程;学会了将实际的问题通过建立模型转化为数学问题（数学方程）并用matlab程序求解结果。

掌握了function函数的使用，并用function函数编辑程序实现二分法迭代和牛顿法迭代输出结果及中间迭代的过程数据，可以说是越用越熟练。训练了数学建模与应用的能力，还训练到了自己梳理程序流程图的能力。

发现了牛顿法与二分法的优缺点。首先，由于二分法在每一步都取区间中点，对于端点取值的大小差别巨大的平滑函数，二分法是极度低效的。其次，对于在区间 [a,b] 中不改变符号的连续函数，如 F(x)=(x−1)2F(x)=(x-1)^2F(x)=(x-1)^2 ，二分法不保证可以求到根。与二分法不同，牛顿法满足二阶收敛，但是牛顿法的二阶收敛率只在局部被满足。特别地是，若初始点离零点不够近，牛顿法可能会发散，或者变得十分低效。虽然对于牛顿法的使用需要注意不收敛的问题，但是对于二次连续可导的凸函数（F″(x)>0，∀x），牛顿法具有很强的收敛性质（单调收敛）。

今后学习中对于题目分析，信息梳理能力有待提高。对于简化编程语言的能力还有特别大的进步空间。

附录：

二分法代码：

f = @(x) x^3+x-3;

format short

[x,i] = bisect(f,1,2,1e-5,1000)

% 输入函数，范围，精度,最大迭代次数

% 输出根,迭代次数

function [x,i] = bisect(f,a,b,eps,max_iter)

    if sign(f(a))~=sign(f(b))

        for i = 1:max_iter  

            c = a/2+b/2;

            disp(i);

            if (b-a)<eps || abs(f(c))<eps

                x = c;

                break

            end

            if sign(f(a))==sign(f(c))

                a = c;

            else

                b = c;

            end

            disp(c);

            disp(f(c));

        end

    end

end

牛顿法代码：

clc

clear all

syms x;

f=input("请输入需求零解的方程f(x)=(自变量为x,如x^3-x^2-5):   ");

p0=input("请输入牛顿迭代法的初始值p_0：  ");

tol=input("请输入精度E：  ");

maxK=input("请输入最大迭代次数：  ");

[p,k,Y]=NTM(f,p0,tol,maxK);

DP=sprintf("使用牛顿迭代法法迭代%d次，计算%s=0以%g为迭代初始值的解为：%g",k,f,p0,p);

disp(DP);

fprintf("迭代值如下：");

disp(Y);

function [p,k,Y]=NTM(f,p0,tol,maxK)

%p0表示迭代初始值

%f表示要求解的方程

%maxK表示规定的最大迭代次数

%tolr表示允许误差

%k表示最终迭代的次数

%p表示最终迭代的值

    syms x;

    P(1)=p0;

    k=2;

    df=diff(f);     %利用diff()函数计算f(x)的导数

    P(k)=P(k-1)-subs(f,x,P(k-1))/subs(df,x,P(k-1)); %第二次迭代的结果

    while k<=maxK

        err=abs(P(k)-P(k-1));    %err表示相邻的迭代值的差值

        if(err<tol)

            fprintf('迭代%d次即可满足允许误差值退出\n',k-1);

            break;

        end

        k=k+1;

        P(k)=P(k-1)-subs(f,x,P(k-1))/subs(df,x,P(k-1));  %迭代

    end      %共迭代了k-1次

    if(k-1==maxK)

        disp("超过最大迭代次数！");

    end

    p=P(k);

    k=k-1;

    Y=P;

end