数值计算方法第2章作业

第二章

一．问题综述：

二．问题分析：

依题意，均质梁的最大挠度值在满足dy

dx

=0时达到，即：

dy dx =

ω0

120EIL

(−5x4+6L2x2−L4)=0

简单进行化简可得需要求解的非线性方程：

−5x4+6L2x2−L4=0

其中，由已知，L=600cm。

观察方程，可得：

−5(x2)2+6L2(x2)−L4=0

这样可以解出精确解，由二次方程求根公式：

x2=

−6L2±√36L4−20L4

解得：x2=360000cm2或72000cm2。由于是一个实际问题，所以显然可知x>0，即可求得x=600cm 或120√5cm≈268.328157299975cm。此即为真值，求解本题时可以利用所求出的真值进行迭代过程中的相对误差的计算。

三．问题解决：

1.二分法

二分法是通过不断取求根区间的中点，然后根据中点的函数值与求根区间两端点的函数值的正负关系来判断根的位置，从而逐渐缩小求根区间，最终求得数值方法下精度允许的解。

源程序bisect.m（定义函数）：

function [c,count,ea,es,yc]=bisect(a,b,tol)

% - a and b are the left and right endpoints

% - tol is the tolerance

% - c is the zero

% - yc=f(c)

% - ea is the err between the calculated value

% - es is the err between the true value

% - count is the times of iteration

f=@(x)(-5*x^4+2*600^2*3*x^2-600^4);

ya=f(a); yb=f(b);

cn=(a+b)/2; %cn is the zero calculated last time

er=1;count=0;

if ya*yb > 0

disp('invaild input');

er=0;

else

while(b-a)>tol %beyond the tolerance

count=count+1;

c=(a+b)/2;

yc=f(c);

if yc==0

a=c;

b=c;

elseif yb*yc>0

b=c;

yb=yc;

else

a=c;

ya=yc;

end

ea=abs((c-cn)/c);

cn=c; %update the zero

es=abs((c-600)/600); %this is the first root

%es=abs((c-268.328157299975)/268.328157299975);%second root fprintf('root=%g ea=%g es=%g\n',c,ea,es);

yc=f(c);

end

c=vpa(c,15);

end

if er==0

c=0;ea=0;es=0;yc=0; %invalid output

end

使用不同的输入值进行迭代，可得：

表1.二分法计算结果

3 590 620 1e-13 47 600.0 3.7896e-16 0 0

4 590 620 1e-10 39 600.000000000018 9.0949e-14 3.0316e-14 -0.0315

5 590 620 1e-5 22 599.999997615814 1.1921e-08 3.9736e-09 4.1199e+03

迭代过程篇幅过长，此处略去，只保留每次使用函数迭代后的输出结果。从迭代过程中可以看出二分法在处理此类两端异号的连续函数时，能够迅速地找到方程根的范围，但是在要求精度较高的情况下时，根的收敛的速度会大大下降，这个现象可以通过资料中的迭代次数公式:

n=⌈log2(x r−x l E a,d

)⌉

求得。但是进行二分法的前提是必须知道每一个根的区间，所以当方程较复杂时，二分法就不是很适用了。从迭代过程中还可以看出，计算值与真实值的相对误差始终小于近似相对误差，这是有相应理论支持的，所以在实际计算过程中，有时可以通过设置近似相对误差限来保证精度要求。

2.试位法

试位法是通过不断求求根区间两端连线与x轴的交点，并比较此交点的函数值与边界点处的函数值，从而判断根的位置，缩小求根区间，最终求得所要求的精度允许的解。

源程序regula.m（定义函数）：

function [c,count,ea,es,yc]=regula(a,b,delta,epsilon)

% - a and b are the left and right endpoints

% - tol is the tolerance of right and left

% - epsilon is the tolerance of yc and 0

% - c is the zero

% - yc=f(c)

% - ea is the err between the calculated value

% - es is the err between the true value

% - count is the times of iteration

f=@(x)(-5*x^4+2*600^2*3*x^2-600^4);

ya=f(a);yb=f(b);

cn=b-yb*(a-b)/(ya-yb); %cn is the zero calculated last time

c=cn;yc=f(c);er=1;count=0;

if ya*yb > 0

disp('invaild input');

er=0;