Columns 13 through 15
1.36878610257799 1.36881787439609 ans =1
1.36880377314363
取解为1.36880377314363
20?2xk2?xk3 对于 xk?1?,只需修改diedai.m文件中的g,把其改为g1,编写m文件
10function y=g1(x); y=(20-2*x^2-x^3)/10;
在命令窗口中输入diedai('g1',1,10^(-9),30) …..
k = 24 er =
1.36601568861328 k = 25 er =
1.36860974051282 k = 26 er =
1.36942327571766
取解为1.36860974051282 牛顿法:编写m文件
function[p1,er,k,y]=ndf(f,df,p0,tol,max) %f是非线性函数 ?是f的微商 %p0是初始值
%tol是给定的允许误差 %max是迭代的最大次数 %p1是牛顿法求得的近似解 %er是p0的误差估计 %k是迭代次数 %y=f(p1) p0,feval('f',p0) for k=1:max
p1=p0-feval('f',p0)/feval('df',p0); er=abs(p1-p0); p0=p1;
p1,er,k,y=feval('f',p1)
if((er 定义函数m文件 function y=f(x) y=x^3+2*x^2+10*x-20; 定义微商函数m文件 function y=df(x) y=3*x^2+4*x+10; 在命令窗口输入 >> ndf('f','df',1,10^(-9),10) p0 = 1 ans = -7 p1 = 1.36933647058824 er = 0.41176470588235 k = 2 y = 0.01114811941245 p1 = 1.36880818861753 K= 3 y= 1.704487072373695e-006 k = 1 y = 0.91756564217382 p1 = er = 0.04242823529412 k = 4 y= 3.907985046680551e-014 p1= 1.36880810782137 er= 1.776356839400251e-015 k= 5 1.41176470588235 p1 = y= 1.36880810782137 0 er = ans= 8.079615732015100e-008 1.36880810782137 由结果知道牛顿法迭代到第三次已经达到所要求的精度 故方程的根为1.36880810782137 3. 编写Steffensen.m文件 i=2;x0=1;%设初始值 f=inline('20/(x^2+2*x+10)');%迭代格式 y=f(x0); z=f(y); x1=x0-(y-x0).^2/(z-2*y+x0); S.result=[x0;x1]; while abs(x1-x0)>=1e-9 %迭代精度 x0=x1; y=f(x0); z=f(y); x1=x0-(y-x0).^2/(z-2*y+x0); i=i+1; S.result(i)=x1; end S.step=[(0:i-1)]'; fprintf('迭代步数为:\\t%d\\n',i-1); for j=1:i fprintf('d',S.step(j));fprintf('\\t'); f 在命令窗口输入Steffensen 迭代步数为: 4 0 1.0000000 1 1.3708139 2 1.3688082 3 1.3688081 4 1.3688081 分析结果知,Steffensen迭代加速步数减少了,到第四步已经达到了精度要求。 4.把迭代格式改为(20-2*x^2-x^3)/10,保存,在命令窗口输入 Steffensen 迭代步数为: 5 0 1.0000000 1 1.3334921 2 1.3684154 3 1.3688081 4 1.3688081 5 1.3688081 分析结果知,Steffensen迭代加速步数减少了,到第五步已经达到了精度要求。 五、用不同的数值方法计算积分 ?3110?10?f(x)dx的近似值,其中f(x)???sin x?x?2⑴ 取不同的步长h,分别用复合梯形公式和复合辛普森公式计算积分,比较两个公式的计 算效果,是否存在一个最小的h,使得精度不能再被改善? ⑵ 用龙贝格求积公式,取??10,并打印出T-表。 解:(1)用复合梯形公式,编写fhtx .m文件 function s=fhtx(f,a,b,n) %f是被积函数 ?分别为积分的上下限 %n是子区间的个数 %s是梯形总面积 h=(b-a)/n; s=0; for k=1:(n-1) x=a+k*h; s=s+feval('f',x); end s=h*(feval('f',a)+feval('f',b))/2+h*s; 编写被积函数文件tf.m function y=f(x); y=(10/x)^2*sin(10/x); 在命令窗口输入 >> fhtx('tf',1,3,10) ans = -4.7789 >> fhtx('tf',1,3,15) ans = -2.8604 >> fhtx('tf',1,3,20) ?4 ans = -2.2208 用复化辛普森公式,编写fhxps.m文件 function s=fhxps(f,a,b,n) %f是被积函数 ?分别为积分的上下限 %n是子区间的个数 %s是梯形总面积 h=(b-a)/(2*n); s1=0;s2=0; for k=1:n x=a+(2*k-1)*h; s1=s1+feval('f',x); end for k=1:n-1 x=a+2*k*h; s2=s2+feval('f',x); end s=h*(feval('f',a)+feval('f',b)+4*s1+2*s2)/3; 在命令窗口输入 >> fhxps('f',1,3,15) ans = -1.4136 >> fhxps('f',1,3,20) ans = -1.4220 取n较大时 >> fhxps('f',1,3,1000) ans = -1.42602475569236 >> fhtx('tf',1,3,1000) ans = -1.42633620719670 >> fhtx('tf',1,3,2000) ans = -1.42610261856845 >> fhxps('f',1,3,2000) ans -1.42602475630535 分析结果知道两种方法有微小误差,当步长越大,精确度就越大,h取2000时两种方法误差精确度到小数点第三位。