1、 數(shù)值計算題目名稱 數(shù)值計算 學 院 機電工程學院 專業(yè)班級 13微電子制造2班 學 號 3113000453 姓 名 肖鎧斌 2016-12-30摘 要 學習數(shù)值計算引論，應用MATLAB編寫Gauss 列主元消元法求解線性方程組、多項式插值、Gauss 積分方法、Euler 方法求常微分方程初值問題、Newton 迭代法求非線性方程的程序。 【關鍵詞】：Gauss 列主元消元法求解線性方程組、多項式插值、Gauss 積分方法、Euler 方法求常微分方程初值問題、Newton 迭代法求非線性方程目錄 1 Gauss 列主元消元法求解線性方程組11.1問題描述11.2計算程序11.3算例22

2、 Lagrange多項式插值32.1問題描述32.2計算程序32.3算例33 Gauss 積分方法43.1.問題描述43.2計算程序43.3算例44 Euler 方法求常微分方程初值問題54.1.問題描述54.2計算程序44.3算例45Newton 迭代法求非線性方程65.1.問題描述65.2計算程序65.3算例7參考文獻8一 、編寫 Gauss 列主元消元法求解線性方程組的程序，要求附有算例。 1、問題描述：計算機中運算的時候常會碰到兩個問題。  1）一旦遇到某個主元等于0，消元過程便無法進行下去。  2）在長期使用中還發(fā)現(xiàn)，即使消元過程能進行下去，但是當某個主元的絕對值

3、很小時，求解出的結果與真實結果相差甚遠。 2、程序： % 輸出的量：系數(shù)矩陣和增廣矩陣的秩RA,RB, 方程組中未知量的個數(shù)n和有關方程組解及其解X的信息.function RA,RB,n,X=gaus(A,b) B=A b; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; %求矩陣A、B的秩，并作比較if zhica>0,disp('RA=RB，此方程組無解.')returnendif RA=RB if RA=n disp('RA=RB=n，此方程組有唯一解.') X=zeros(n,1); C=zero

4、s(1,n+1); %生成全零矩陣 for p= 1:n-1 Y,j=max(abs(B(p:n,p); %返回行向量Y和j，y向量記錄abs(B(p:n,p) 的每列的最大值，j向量記錄每列最大值的行號C=B(p,:); %把矩陣B中第p行所有列的值全都賦給矩陣CB(p,:)= B(j+p-1,:); B(j+p-1,:)=C; %交換A,b第jk與第k 行元素 for k=p+1:n m= B(k,p)/B(p,p); %消元計算 B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1); endend b=B(1:n,n+1); A=B(1:n,1:n); X(n)=b(

5、n)/A(n,n); for q=n-1:-1:1 %從n-1循環(huán)到1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); %回代計算 endelse disp('RA=RB<n，此方程組有無窮多解.') endEnd 3、算例：求解線性方程組在命令窗口輸入：>> A=5 2 1;2 8 -3;1 -3 -6; b=8;21;1; RA,RB,n,X =gaus (A,b)RA=RB=n，此方程組有唯一解.RA = 3RB = 3n = 3X = 1 2-1二 編寫多項式插值的程序，要求附有算例 1、問題描述：Lagrange

6、差值法滿足插值條件的、次數(shù)不超過n的多項式是存在而且是唯一的。 2、程序：function yy=lagrange(x,y,x0) /定義拉格朗日差值函數(shù) xlen=length(x); /X向量的維數(shù) ylen=length(y); /y向量的維數(shù)if xlen=ylen warning('The length of x and y is not equal'); /提示x.y維數(shù)不等endn=length(x0); for i=1:n /循環(huán)。從1到n z=x0(i); s=0.0; for k=1:xlen p=1.0; for j=1:ylen if j=k p=p*(

7、z-x(j)/(x(k)-x(j); end end s=p*y(k)+s; end yy(i)=s; end 3、算例：已知數(shù)據(jù)如下表，試用lagrange差值多項式求x分別為0.5626;0.5635;0.5645時函數(shù)的近似值。0.561600.562800.564010.565210.827410.826590.825770.81495解：在命令窗口輸入執(zhí)行命令，可得結果：>> x=0.56160;0.56280;0.56401;0.56521;y=0.82741;0.82659;0.82577;0.81495;x0=0.5626;0.5635;0.5645;y0=lagr

8、ange(x,y,x0)y0 = 0.8265 0.8268 0.8231>> plot(x,y,'o',x0,y0,'k')得到圖形三 編寫 Gauss 積分方法的程序，要求附有算例。 1、問題描述：被積分函數(shù)f(，)一般是很復雜的，即使能夠得出它的顯式，其積分也是很繁的。因此，一般用數(shù)值積分來代替函數(shù)的定積分。 2、程序：% f：積分表達式，可以是函數(shù)句柄、inline函數(shù)、匿名函數(shù)、字符串表達式，但是必須可以接受矢量輸入% a,b：積分上下限，注意積分區(qū)間太大會降低精度，此時建議使用復化求積公式，默認-1 1% n：積分階數(shù)，可以任意正整數(shù)，但

9、是不建議設置過大，大不一定能得到更好的精度，默認7% tol：積分精度，默認1e-6% ql：積分結果% Ak：求積系數(shù)% xk：求積節(jié)點，滿足ql=sum(Ak.*fun(xk)function ql,Ak,xk=guass(f,a,b,n,tol)if nargin=1 % 輸入?yún)?shù)有效性檢驗 a=-1;b=1;n=7;tol=1e-8;elseif nargin=3 n=7;tol=1e-8;elseif nargin=4 tol=1e-8;elseif nargin=2|nargin>5 error('The Number of Input Arguments Is Wr

10、ong!');endsyms x % 計算求積節(jié)點p=sym2poly(diff(x2-1)(n+1),n+1)/(2n*factorial(n);tk=roots(p); % 求積節(jié)點Ak=zeros(n+1,1); % 計算求積系數(shù)for i=1:n+1 xkt=tk; xkt(i)=; pn=poly(xkt); fp=(x)polyval(pn,x)/polyval(pn,tk(i); Ak(i)=quadl(fp,-1,1,tol); % 求積系數(shù)endxk=(b-a)/2*tk+(b+a)/2; % 積分變量代換，將a,b變換到-1,1f=fcnchk(f,'vec

11、torize'); % 檢驗積分函數(shù)fun有效性fx=f(xk)*(b-a)/2; % 計算變量代換之后積分函數(shù)的值ql=sum(Ak.*fx); % 計算積分值 3、算例，求的值。編寫函數(shù)文件：Gaussf.mfunction f=gaussf(x)f=cos(x);End在命令窗口輸入：>> ql,Ak,xk=guass('gaussf', 0, 1, 7, 1e-6)ql = 0.8415Ak = 0.1012 0.2224 0.1012 0.2224 0.3137 0.3137 0.3627 0.3627xk = 0.0199 0.1017 0.98

12、01 0.8983 0.2372 0.7628 0.4083 0.5917四 編寫 Euler 方法求常微分方程初值問題的程序，要求附有算例。 1、問題描述：歐拉法的特點：單步，顯式，一階求導精度，截斷誤差貳階 2、程序：歐拉公式：Euler.mfunction xx,yy=euler(f,a,b,y0,h) x=a:h:bn=length(x);y=zeros(1,n);y(1)=y0;for i=1:n-1; y(i+1)=y(i)+h*feval(f,x(i),y(i);end xx=x; yy=y;plot(xx,yy,'rp')xlabel('Independ

13、ent Variable x'),ylabel('Dependent Variable y')title('Solution of ODE with modified euler method') 3、算例：用歐拉方法解初值問題 ，取h=0.1.解：建立函數(shù)文件myfun.mfunction f=myfun(x,y)f=1+x2+y2End在窗口輸入命令：>>x y=euler('myfun',0,1,0,0.1)x = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.

14、8000 0.9000 1.0000y = 0 0.1000 0.2020 0.3101 0.4287 0.5631 0.7198 0.9076 1.1390 1.4327 1.8189所得圖形： 改進后的歐拉公式：function xx,yy=euler(f,a,b,y0,h)x=a:h:bn=length(x);y=zeros(1,n);y(1)=y0;for i=1:n-1; yp=y(i)+h*feval(f,x(i),y(i); yc=y(i)+h*feval(f,x(i+1),yp); y(i+1)=(yp+yc)/2; end xx=x; yy=y;ey=-1./xplot(xx

15、,yy,'rp')xlabel('Independent Variable x'),ylabel('Dependent Variable y')title('Solution of ODE with modified euler method')在窗口輸入命令：>>x y=euler('myfun',0,1,0,0.1)x = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000y = 0 0.1010 0.206

16、1 0.3196 0.4469 0.5943 0.7710 0.9899 1.2719 1.6529 2.2020五 編寫 Newton 迭代法求非線性方程的程序，要求附有算例。 1、問題描述：構造迭代函數(shù)的一條重要途徑是用近似方程來代替原方程去求根。因此，如果能將非線性方程f x =0用線性方程去代替，那么，求近似根問題就很容易解決，而且十分方便。牛頓(Newton)法就是一種將非線性方程線化的一種方法 2、程序：Newton.m%y=newton(a,n,x0,nn,eps1) 輸入變量： %a n+1元素的一維實數(shù)組，輸入?yún)?shù)，按升冪存放

17、方程系數(shù)。 %n 整變量，輸入?yún)?shù)，方程階數(shù)。 %x0 實變量，輸入?yún)?shù)，初始迭代值。 %nn 整變量，輸入?yún)?shù)，允許的最大迭代次數(shù)。 %eps1 實變量，輸入?yún)?shù)，控制根的精度。function y=newton(a,n,x0,nn,eps1)x(1)=x0; %初始向量b=1;i=1;while(abs(b)>eps1*x(i) %abs(b)取模，判斷 i=i+1;x(i)=x(i-1)-n_f(a,n,x(i-1)/n_df(a,n,x(i-1); %b=x(i)-x(i-1);if(i>nn) error('nn is full'); return;endendy=x(i);i程序中調兩個子函數(shù)n_f.m和n_df.m文件如下： n_f.m: function y=n_f(a,n,x)%待求根的實數(shù)代數(shù)方程的函數(shù) y=0.0; for i=1:(n+1) y=y+a(i)*x(n+1-i); end  n_df.m: