数值分析实验报告三
- 格式：doc
- 大小：169.50 KB
- 文档页数：15
实验目的 和 要 求
实验内容和步骤：
实验内容：
1、Compare the number of computation for finding the root of
e x 10 x 2 0 with accuracy 105 .
(1)、Use the bisection method starting with the interval [0,1]. (2) 、 Use the iteration method xk 1 (2 ex ) /10 ,the initial value
数学与信息工程学院
实 验 报 告
课程名称： 实 验 室： 实验台号： 班 姓 级： 名：
计算方法 7404
实验日期：
2014
年 5
月 21 日
实验名称
非线性方程求解及方程组直接解法 1.二分法的 Matlab 实现； 2.牛顿法的 Matlab 实现； 3.牛顿下山法、割线法、艾特金加速法、重根迭代法、非线性方程 组牛顿法中任选其ห้องสมุดไป่ตู้。 4.列主元消元法的 Matlab 实现； 5.二选一： （1）平方根法的 Matlab 实现； （2）追赶法的 Matlab 实 现。
2. 建立名为fun.m的M文件
function fun=fun(x) fun=exp(x)+10*x-2;
3.在 MATLAB 工作窗口输入：
x=-4:0.1:4; y=exp(x) +10*x-20; plot(x,y) grid [k,x,wuca,yx]=erfen (﹣1,1,10^-5)
实验步骤： 1） 2） 实验编程 运行结果 二分法 1） 实验编程
1. 新建 M 文件 newtonqx.m 输入： function [k,x,wuca,yx]=erfen(a,b,abtol) a(1)=a; b(1)=b; ya=fun(a(1)); yb=fun(b(1)); %程序中调用的fun.m 为函数 if ya* yb>0, disp('注意：ya*yb>0,请重新调整区间端点a和b.'), return end max1=-1+ceil((log(b-a)- log(abtol))/ log(2)); % ceil(u)是大于u的最小 取整数 for k=1: max1+1 a;ya=fun(a); b;yb=fun(b); x=(a+b)/2; yx=fun(x); wuca=abs(b-a)/2; k=k-1; [k,a,b,x,wuca,ya,yb,yx] if yx==0 a=x; b=x; elseif yb*yx>0 b=x;yb=yx; else a=x; ya=yx; end if b-a< abtol , return, end end k=max1; x; wuca; yx=fun(x);
2）运行结果
ans = 0 ans = 1.0000 ans = 2.0000 ans = 3.0000 ans = 4.0000 ans = 5.0000 ans = -1.0000 0 0 0 0 0.0625 1.0000 1.0000 0.5000 0.2500 0.1250 0.1250 0 0.5000 0.2500 0.1250 0.0625 0.0938 1.0000 0.5000 0.2500 0.1250 0.0625 0.0313 -11.6321 -1.0000 -1.0000 -1.0000 -1.0000 -0.3105 10.7183 10.7183 4.6487 1.7840 0.3831 0.3831 -1.0000 4.6487 1.7840 0.3831 -0.3105 0.0358
-0.3105 -0.1375 -0.0509 -0.0076 -0.0076 -0.0076 -0.0021 -0.0021 -0.0008 -0.0001 -0.0001 -0.0001
0.0358 0.0358 0.0358 0.0358 0.0141 0.0033 0.0033 0.0006 0.0006 0.0006 0.0002 0.0001
3、用迭代法求方程 x3 x2 1 0在x0 1.5 附近的一个根。方程写成下 列等价形式，并建立相应的迭代公式。
x 1 1 1 ，迭代公式 xk 1 1 2 2 x xk
4、用列主元消法解线性方程组 Ax=b。
3 2 6 4 A 10 7 0 , b 7 5 1 5 6
x(i+1)=x(i)-fnq(x(i))/(dfnq(x(i))+eps); piancha=abs(x(i+1)-x(i)); xdpiancha= piancha/( abs(x(i+1))+eps); i=i+1; xk=x(i);yk=fnq(x(i)); [(i-1) xk yk piancha xdpiancha] if (abs(yk)<ftol)&((piancha<tol)|(xdpiancha< tol)) k=i-1; xk=x(i);[(i-1) xk yk piancha xdpiancha] return; end end if i>gxmax disp('请注意：迭代次数超过给定的最大值gxmax。') k=i-1; xk=x(i);[(i-1) xk yk piancha xdpiancha] return; end [(i-1),xk,yk,piancha,xdpiancha]'; 2. 建立名为 fnq.m 的 M 文件 function y=fnq(x) y=2*x.^4+24*x.^3+61*x.^2-16*x+1; 3. 建立名为 dfnq.m 的 M 文件 function y=dfnq(x) y=8*x.^3+72*x.^2+122*x-16; 4. 在 MATLAB 工作窗口输入： >>x=0.5:0.1:4; >> y=2*x.^4+24*x.^3+61*x.^2-16*x+1; >> plot(x,y) >>grid >> [k,xk,yk,piancha,xdpiancha]=newtonqx(0.5,0.001, 0.001,100)
3.0037 0.7752 0.1972 0.0497 0.0125 0.0031 0.0008 0.0002 0.0000 0.0000
0.1777 0.0967 0.0508 0.0260 0.0132 0.0066 0.0033 0.0016 0.0007 0.0007
0.5515 0.4287 0.2903 0.1751 0.0973 0.0514 0.0262 0.0129 0.0056 0.0056
6.0000 ans = 7.0000 ans = 8.0000 ans = 9.0000 ans = 10.0000 ans = 11.0000 ans = 12.0000 ans = 13.0000 ans = 14.0000 ans = 15.0000 ans = 16.0000 ans = 17.0000
-0.1375 -0.0509 -0.0076 0.0141 0.0033 -0.0021 0.0006 -0.0008 -0.0001 0.0002 0.0001 -0.0000
k =17 x =0.0905 wuca =7.6294e-006 yx =-2.5908e-005
牛顿法 1）实验编程
6. 在 MATLAB 工作窗口输入： x0=1.5; k=3; diedai1(x0,k)
0.0625 0.0781 0.0859 0.0898 0.0898 0.0898 0.0903 0.0903 0.0905 0.0905 0.0905 0.0905
0.0938 0.0938 0.0938 0.0938 0.0918 0.0908 0.0908 0.0906 0.0906 0.0906 0.0905 0.0905
2）
运行结果
ans = 1.0000 0.3223 ans = 2.0000 0.2256 ans = 3.0000 0.1748 ans = 4.0000 0.1488 ans = 5.0000 0.1356 ans = 6.0000 0.1289 ans = 7.0000 0.1256 ans = 8.0000 0.1240 ans = 9.0000 0.1233 ans = 9.0000 0.1233 k=9 xk = 0.1233 yk = 3.4012e-005
0.0781 0.0859 0.0898 0.0918 0.0908 0.0903 0.0906 0.0905 0.0905 0.0905 0.0905 0.0905
0.0156 0.0078 0.0039 0.0020 0.0010 0.0005 0.0002 0.0001 0.0001 0.0000 0.0000 0.0000
1. 新建 M 文件 newtonqx.m 输入：
%输入：初始值x0;近似值xk 的精度tol； f(xk)的精度tol。 function [k,xk,yk,piancha,xdpiancha]=newtonqx(x0,tol,ftol,gxmax) x(1)=x0; for i=1: gxmax
piancha = 6.9657e-004 xdpiancha = 0.0056
埃特金加速收敛方法
1）实验编程 1. 新建 M 文件 diedail.m 输入：
function [k,piancha,xdpiancha,xk]=diedai1(x0,k) % 输入的量--x0是初始值,k是迭代次数 x(1)=x0; for i=1:k x(i+1)=fun1(x(i));%程序中调用的fun1.m为函数y=φ (x) piancha= abs(x(i+1)-x(i)); xdpiancha=piancha/( abs(x(i+1))+eps); i=i+1;xk=x(i);[(i-1) piancha xdpiancha xk] end if (piancha >1)&(xdpiancha>0.5)&(k>3) disp('注意：此迭代序列发散，请重新输入新的迭代公式') return end if (piancha < 0.001)&(xdpiancha< 0.0000005)&(k>3) disp('祝贺您！此迭代序列收敛，且收敛速度较快') return end p=[(i-1) piancha xdpiancha xk]'; 2. 建立名为 fnq.m 的 M 文件 function f=fnq(x) f=x.^3+2*x.^2+10*x-20; 3. 建立名为 dfnq.m 的 M 文件 function f=dfnq(x); f=3*x^2+4*x+10; 4. 建立名为 fai.m 的 M 文件 function f=fai(x); f=x-fnq(x)/(dfnq(x)+eps); 5. 建立名为 fun1.m 的 M 文件 function f=fun1(x) f=x-(fai(x)-x)^2/(fai(fai(x))-2*fai(x)+x);
