二分法非线性方程求解.docx
《二分法非线性方程求解.docx》由会员分享，可在线阅读，更多相关《二分法非线性方程求解.docx（17页珍藏版）》请在冰点文库上搜索。
二分法非线性方程求解
1、编程实现以下科学计算算法，并举一例应用之（参考书籍《精通MATLAB科学计算》，王正林等编著，电子工业出版社，2009年）
“二分法非线性方程求解”
二分法的具体求解步骤如下。
（1）计算函数f（x）在区间[a,b]中点的函数值f（（a+b）/2），并作下面的判断：
如果
，转到
（2）；
如果
，令
，转到
（1）；
如果
，则
为一个跟。
（2）如果
（
为预先给定的精度），则
为一个根，否则令
，转到
（1）。
在MATLAB中编程实现的二分法函数为：
HalfInterval。
功能：
用二分法求函数在某个区间上的一个零点。
调用格式：
root=HalfInterval（f,a,b,eps）.
其中，f函数名；
a为区间左端点；
b为区间右端点；
eps为根的精度；
root为求出的函数零点。
二分法的MATLAB程序代码如下：
functionroot=HalfInterval（f,a,b,eps）
%二分法求函数f在区间[a,b]上的一个零点
%函数名：
f
%区间左端点：
a
%区间右端点：
b
%根的精度：
eps
%求出的函数零点：
root
if（nargin==3）
eps=1.0e-4;
end
f1=subs（sym（f）,findsym（sym（f））,a）;%两端点的函数值
f2=subs（sym（f）,findsym（sym（f））,b）;
if（f1==0）
root=a;
end
if（f2==0）
root=b;
end
if（f1*f2>0）
disp（'两端点函数值乘积大于0!
'）;
return;
else
root=FindRoots（f,a,b,eps）;%调用求解子程序
end
functionr=FindRoots（f,a,b,eps）
f_1=subs（sym（f）,findsym（sym（f））,a）;
f_2=subs（sym（f）,findsym（sym（f））,b）;
mf=subs（sym（f）,findsym（sym（f））,（a+b）/2）;%中点函数值
if（f-1*mf>0）
t=（a+b）/2;
r=FindRoots（f,t,b,eps）;%右递归
else
if（f_1*mf==o）
r=（a+b）/2;
else
if（abs（b-a）<=eps）
r=（b+3*a）/4;%输出根
else
s=（a+b）/2;
r=FindRooots（f,a,b,eps）;%左递归
end
end
end
流程图：
实例应用：
采用二分法求方程
在区间[0，1]上的一个根。
解：
流程图：
在MATLAB命令窗口中输入：
r=HalfInterval（'x^3-3*x+1',0,1）
运行结果：
2、编程以解决以下科学计算问题。
试验6现有一平面上的封闭曲线，取一点建立坐标系，每隔
弧度测一点，数据如下表：
i
0和18
1
2
3
4
5
6
7
8
Xi
100
134
164
180
198
195
186
160
136
Yi
503
525
514.3
451.0
326.5
188.6
92.2
59.6
62.2
i
9
10
11
12
13
14
15
16
17
Xi
100
66
35
15
0
5
17
32
63
yi
102.7
147.1
191.6
236.0
280.5
324.9
324.9
413.8
458.3
用周期样条求曲线轮廓并作图。
分析：
将周期样条分成两部分来求，最后画在一个图上。
用spline进行拟合。
流程图：
源程序：
>>x1=[05173263100134164180198];
y1=[280.5324.9369.4413.8458.3503525514.3451326.5];
x11=[0:
1:
198];
y11=spline（x1,y1,x11）;
plot（x11,y11,'*-',x1,y1,'-.rd'）
holdon
x2=[1981951861601361006635150];
y2=[326.5188.692.259.662.2102.7147.1191.6236.0280.5];
x22=[198:
-1:
0];
y22=spline（x2,y2,x22）;
plot（x22,y22,'*-',x2,y2,'-.rd'）
legend（'计算数据','实验数据'）
运行结果：
实验7用电压V=10伏的电池给电容器，电容器上t时刻的电压
，其中
是电容器的初始电压，
是充电常数，试由下面一组
，
数确定
和
。
t/s
0.5
1
2
3
4
5
6
7
V/V
6.36
6.48
7.26
8.22
8.66
8.99
9.43
9.63
分析：
采用最小二乘法进行拟合，对V（t）=V-（V-V0）exp（-t/τ）两边求自然对数得到：
log（V-V（t））=log（V-V0）-t/τ
令k=-1/τ，w=log（V-V0），x=t，y=log（V-V（t））。
得到方程：
y=kx+w
流程图：
编程如下：
t=[0.51234579];
V=[6.366.487.268.228.668.999.439.63];
V1=log（10-V）;
[p,s]=polyfit（t,V1,1）;
[V1p,delta]=polyval（p,t,s）;
t0=-1/p
（1）
V0=10-exp（p
（2））
plot（t,V,'--r',t,10-（10-V0）*exp（-t/t0））
legend（'实验数据','拟合数据'）
结果：
t0=
3.5269
V0=
5.6221
实验8假定某天的气温变化记录如下表，试用最小二乘方法找出这一天的气温变化规律。
t/h
0
1
2
3
4
5
6
7
8
9
10
11
12
T/
15
14
14
14
14
15
16
18
20
22
23
25
28
t/h
13
14
15
16
17
18
19
20
21
22
23
24
T/
31
32
31
29
27
25
24
22
20
18
17
16
考虑下列类型函数，计算误差平方和，并作图比较效果：
（1）二次函数；
（2）三次函数；
（3）四次函数；
（4）函数
。
分析：
对“二次函数”“三次函数”“四次函数”用最小二乘法进行拟合（polyfit）,对指数函数首先两边取对数，再通过移项化简，获得普通形式的函数，用最小二乘法进行拟合，并求出相应系数。
流程图：
源程序
t=0:
1:
24;
T=[15141414141516182022232528313231292725242220181716];
%二次函数拟合
[p2,s2]=polyfit（t,T,2）;
T2=polyval（p2,t）;
plot（t,T,'*-',t,T2）;
legend（'观测数据','计算数据'）
title（'二次函数拟合'）
p2
deltaT2=sum（（T2-T）.*（T2-T））
%三次函数拟合
[p3,s3]=polyfit（t,T,3）;
T3=polyval（p3,t）;
figure
plot（t,T,'*-',t,T3）;
legend（'观测数据','计算数据'）
title（'三次函数拟合'）
p3
deltaT3=sum（（T3-T）.*（T3-T））
%四次函数拟合
[p4,s4]=polyfit（t,T,4）;
T4=polyval（p4,t）;
figure
plot（t,T,'*-',t,T4）;
legend（'观测数据','计算数据'）
title（'四次函数拟合'）
p4
deltaT4=sum（（T4-T）.*（T4-T））
%指数函数拟合
Te0=log（T）;
[pe,se]=polyfit（t,Te0,2）;
b=pe
（1）;
c=pe
（2）/2/b;
a=exp（pe（3）+c）;
Te1=polyval（pe,t）;
Te2=exp（Te1）;
figure
plot（t,T,'*-',t,Te2）;
legend（'观测数据','计算数据'）
title（'指数函数拟合'）
a
b
c
deltaTe=sum（（Te2-T）.*（Te2-T））
结果：
二次函数拟合：
三次函数拟合：
四次函数拟合：
指数函数拟合：
误差平方和：
二次方
三次方
四次方
幂函数
误差平方和
241.2430
106.0776
36.2838
178.6060
因此，可以看出，四次函数拟合的最好，其次是三次函数，二次函数拟合的最差。
升级会员 