數(shù)值分析試題

1、1.設(shè)有某實驗數(shù)據(jù)如下：0.10.20.30.40.50.60.70.80.95.12345.30535.56845.93786.42707.07987.94939.025310.3627（1）在MATLAB中作圖觀察離散點(diǎn)的結(jié)構(gòu)，用最小二乘法擬合一個合適的多項式函數(shù)；（2）在MATLAB中作出擬合曲線圖. 解：（1）在MATLAB中作圖觀察離散點(diǎn)的結(jié)構(gòu)，用最小二乘法擬合一個合適的多項式函數(shù)。具體程序如下：>> x=0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9;>> y=5.1234 5.3053 5.5684 5.9378 6.4270 7.

2、0798 7.9493 9.0253 10.3627;>> plot(x,y,'r+'),legend('x,y'),xlabel('x'),ylabel('y'),title('(x,y)的散點(diǎn)結(jié)構(gòu)圖')從圖中各點(diǎn)可以看到，（x,y）散點(diǎn)結(jié)構(gòu)圖的變化趨勢與二次多項式很接近，故選取曲線函數(shù)為二次多項式函數(shù)，擬合多項式的次數(shù)為2。用最小二乘法擬合一個合適的多項式函數(shù)，具體程序如下：>> s=polyfit(x,y,2);>> P=poly2str(s,'t')P =

3、 8.2191 t2 - 1.8822 t + 5.3139則曲線函數(shù)為發(fā)f（x）=8.2192t2-1.8822t+5.3139（2） 在MATLAB中作出擬合曲線圖，具體程序如下：>> x=linspace(0.1,0.9,9);>> y=8.2191.*x.2-1.8822.*x+5.3139;>> plot(x,y),legend('x,y'),xlabel('x'),ylabel('y'),title('擬合曲線y=f(x)')2. 在MATLAB中用復(fù)合Simpson公式編程計算下列

4、積分，并使其結(jié)果至少有6位有效數(shù)字. 解：（1）所安裝的MATLAB中不含有復(fù)合辛普森程序，用M文件編輯器編寫被積函數(shù)的M文件，并以jfSimpson.m命名保存至MATLAB搜索路徑下。程序如下：function I,step = jfSimpson(f,a,b,type,eps)%type = 1 辛普森公式%type = 2 復(fù)合辛普森公式 if(type=2 && nargin=8) eps=1.0e-8; %缺省精度為0.0001end I=0;switch type case 1, I=(b-a)/6)*(subs(sym(f),findsym(sym(f),a)+

5、. 4*subs(sym(f),findsym(sym(f),(a+b)/2)+. subs(sym(f),findsym(sym(f),b); step=1; case 2, n=2; h=(b-a)/2; I1=0; I2=(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)/h; while abs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/6)*(subs(sym(f),findsym(sym(f

6、),x)+. 4*subs(sym(f),findsym(sym(f),(x+x1)/2)+. subs(sym(f),findsym(sym(f),x1); end end I=I2; step=n;end （2）在命令窗口中，調(diào)用jfSimpson函數(shù)進(jìn)行求解。命令窗口顯示如下：>> I,s=jfSimpson(inline('exp(x.2)'),0,2,2,1.0e-8)I = 16.4526s = 102說明：求解積分函數(shù)為16.4526。3. 在MATLAB中試用經(jīng)典四階Runge-Kutta法編程求解初值問題 解：（1）用M文件編輯器將RK4.m、RK

7、_fun.m、RK_main.m的主程序命名并保存至MATLAB搜索路徑下。主程序如下：RK4.mfunction t,y = RK4(func,t0,tt,y0,N,varagin)% Rk方法計算一階級微分方程組，% 由微分方程知識可以知道，對于高階微分方程可以化為一階微分方程組。% 初始時刻為t0，結(jié)束時為tt，初始時刻函數(shù)值為y0% N為 步數(shù) if nargin<4 N=100 end % 如果沒有輸入N的話，那么默認(rèn)計算步長為（tt-t0）/100 y(1,:) = y0(:)'h = (tt-t0)/N; %步長t = t0+0:N'*h; for k =

8、1:N f1 = h*feval(func,t(k),y(k,:); f1 = f1(:)'% RK中的k1f2 = h*feval(func,t(k) + h/2,y(k,:) + f1/2); f2 = f2(:)'% RK中的k2 f3 = h*feval(func,t(k) + h/2,y(k,:) + f2/2); f3 = f3(:)'% RK中的k3f4 = h*feval(func,t(k) + h,y(k,:) + f3); f4 = f4(:)'% RK中的k4 y(k + 1,:) = y(k,:) + (f1 + 2*(f2 + f3)

9、+ f4)/6; % 所求結(jié)果endRK_fun.m%RK4方法的測試函數(shù)function f=RK_fun(t,y)f=-3*y2+2*t2+3*t;RK_main.m%rk方法的主函數(shù) x0=0;xt=1;y0=1;%符號計算給出分析解y=dsolve('Dy=-3*y2+2*t2+3*t','y(0)=1','t') %RK4方法給出計算結(jié)果x,yrk = RK4('RK_fun',x0,xt,y0,10); %對應(yīng)于真解的離散數(shù)據(jù) yreal=subs(y,x); tol=yreal-yrk; %RK4方法與分析解的誤差p

10、lot(x,yreal,x,yrk,'+',x,tol,'*')legend('分析解','RK4計算結(jié)果','兩者誤差') yrk;x,yrk,tol%yrk為所要計算的結(jié)果（2） 運(yùn)行結(jié)果如下： >> RK_main y =<1x1 sym>ans= <1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym>

11、<1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym>

12、<1x1 sym><1x1 sym><1x1 sym><1x1 sym><1x1 sym>4. 在MATLAB中利用Newton法編程計算方程 的一個根. 解：(1)所安裝的MATLAB中不含有newton.m程序，用M文件編輯器編寫被積函數(shù)的M文件，并以newtonqxf.m命名保存至MATLAB搜索路徑下。程序如下：function root=newtonqxf(f,a,b,eps)% 牛頓法求函f數(shù)在區(qū)間a,b上的一個零點(diǎn)% f 為函數(shù)名% a 為區(qū)間左端點(diǎn)% b 為區(qū)間右端點(diǎn)% eps 為根的精度% root 為求出的函數(shù)零點(diǎn)

13、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;endif(f2=0) root=b;end if(f1*f2>0) disp('兩端點(diǎn)函數(shù)值乘積大于0!'); return;else tol=1; fun=diff(sym(f); fa=subs(sym(f),findsym(sym(f),a); fb=subs(sym(f),findsym(sym(f),b); dfa=subs(sym(fun),

14、findsym(sym(fun),a); dfb=subs(sym(fun),findsym(sym(fun),b); if(dfa>dfb) root=a-fa/dfa; else root=b-fb/dfb; end while(tol>eps) r1=root; fx=subs(sym(f),findsym(sym(f),r1); dfx=subs(sym(fun),findsym(sym(fun),r1); root=r1-fx/dfx; tol=abs(root-r1); endend (2) 在命令窗口中，調(diào)用newtonqxf函數(shù)進(jìn)行求解。命令窗口顯示如下:>&

15、gt; root=newtonqxf('x-log(x)-2',3,1.0e-9)root =3.14625. 已知線性方程組 （1）利用MATLAB說明解該方程組的Gauss-Seidel迭代法是否收斂； （2）在MATLAB中利用Gauss-Seidel迭代法求解該方程組，并給出迭代次數(shù). 解：通過M文件編輯器加入gauseidel.m函數(shù)程序，如下：function x,n=gauseidel(A,b,x0,eps,M)% 求解線性方程組的迭代法,其中,% A為方程組的系數(shù)矩陣;% b為方程組的右端項;% x0為迭代初始化向量% eps為精度要求,缺省值為1e-5;% M

16、為最大迭代次數(shù),缺省值100;% x為方程組的解;% n為迭代次數(shù);if nargin=3 eps= 1.0e-6; M = 200;elseif nargin = 4 M = 200;elseif nargin<3 error return;end D=diag(diag(A); %求A的對角矩陣L=-tril(A,-1); %求A的下三角陣U=-triu(A,1); %求A的上三角陣G=(D-L)U;f=(D-L)b;x=G*x0+f;n=1; %迭代次數(shù)while norm(x-x0)>=eps x0=x; x=G*x0+f; n=n+1; if(n>=M) disp('Warning: 迭代次數(shù)太多，可能不收斂！'); return; endend（1）在命令窗口中，輸入程序如下：>> A=5,2,1;-1,4,2;2,-5,10A = 5 2 1 -1 4 2 2 -5 10>> D=diag(diag(A)D = 5 0 0 0 4 0 0 0 10>> L=-tril(A,1)L = -5 -2 0 1 -4 -2 -2 5 -10>> U=-triu(A,1)U = 0 -2 -1 0 0 -2 0 0 0>> G=(D-L)UG = 0 -0.1945 -