數(shù)值計算方法實驗報告

數(shù)值計算方法實驗報告葉耀北師珠2016.01.10&1：實驗一:基本變量與運算：for循環(huán)：s=0fori=1:n（n為常數(shù)）s=s+i;ends函數(shù)及函數(shù)實驗方法：分段函數(shù)的例子：function[y]=func(x)ifx>0y=sin(x);elsey=1-cos(x);end基本繪圖命令實驗：下面看如何畫出這個分段函數(shù)的圖形：x=-pi:0.1:pi;y=[];fori=1:length(x)y(i)=func(x(i));endplot(x,y)畫圖的程序：x=-5:0.1:5;y=[];fori=1:length(x)y(i)=func2(x(i));endplot(x,y)gridon&2：lagrange插值多項式實現(xiàn)：完善的拉格朗日插值基函數(shù)的程序：function[f]=language(x,y,x0)%theLanguageploynomialfunction%thevectorxandyareknownsymst;if(length(x)==length(y))n=length(x);elsedisp('xandyhavediffierentdimensions!')return;endf=0;for(i=1:n)l=y(i);for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));endfor(j=i+1:n)l=l*(t-x(j))/(x(i)-x(j));endf=f+l;simplify(f);if(i==n)if(nargin==3)f=subs(f,'t',x0);elsef=collect(f);f=vpa(f,6);endendendend&3：Newton插值多項式的程序：function[y,R,A,C,L]=newdscg(X,Y,x,M)n=length(X);m=length(x);fort=1:mz=x(t);A=zeros(n,n);A(:,1)=Y';s=0.0;p=1.0;q1=1.0;c1=1.0;forj=2:nfori=j:nA(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1));endq1=abs(q1*(z-X(j-1)));c1=c1*j;endC=A(n,n);q1=abs(q1*(z-X(n)));fork=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k);endy(k)=polyval(C,z);endR=M*q1/c1;L(k,:)=poly2sym(C)&4：龍貝格積分法的上機實現(xiàn)：formatlongf=input('請輸入原函數(shù)f=','s');a=input('積分下限a=');b=input('積分上限b=');eps1=input('精度eps1=');T(1)=double((b-a)/2*(limit(sym(f),findsym(sym(f)),a)+limit(sym(f),findsym(sym(f)),b)));fork=2:4sum1=0;fori=1:2^(k-2)sum1=sum1+subs(sym(f),findsym(sym(f)),(a+(2*i-1)*(b-a)/2^(k-1)));end

T(k)=1/2*T(k-1)+(b-a)/(2^(k-1))*sum1;endfork=1:3S(k)=T(k+1)+1/(4-1)*(T(k+1)-T(k));endfork=1:2c(k)=S(k+1)+1/(4^2-1)*(S(k+1)-S(k));endR(1)=C(2)+1/(4^3-1)*(C(2)-C(1));k=3;while1T(1)=T(2);T(2)=T(3);T(3)=T(4);sum2=0;fori=1:2^ksum2=sum2+subs(sym(f),findsym(sym(f)),(a+(2*i-1)*(b-a)/2^(k+1)));end

T(4)=1/2*T(4)+(b-a)/2^(k+1)*sum2;S(1)=S(2);S(2)=S(3);S(3)=T(4)+1/(4-1)*(T(4)-T(3));C(1)=C(2);C(2)=S(3)+1/(4^2-1)*(S(3)-S(2));R(2)=C(2)+1/(4^3-1)*(C(2)-C(1));ifabs(R(2)-R(1))<eps1break;endR(1)=R(2);endvpa(R(2),9)&5：總結(jié)matlab自己帶有的插值、積分、解線性方程組、解非線性方程及方程組的命令，用實例來展示各命令的用法:&1：Matlab基礎(chǔ)應(yīng)用的實例：就給出美國1930年到2000年的人口數(shù)，根據(jù)已有的數(shù)據(jù)，尋找一種辦法來預(yù)測一下2000年的美國人口數(shù)：19301940195019601970198019901.232031.316691.506791.793232.032122.265052.49633解：clcx=1940:10:1990;y=[1.232031.316691.506791.793232.032122.265052.49633];x0=2000;y0=lagrange(x,y,x0)結(jié)果為：y0=3.0394&2Lagrange插值多項式的程序?qū)嵗簒=1:4;y=x.^2;lagrange(x,y,1)ans=1lagrange(x,y,2)ans=4lagrange(x,y,3)ans=9&3Newton插值多項式的程序?qū)嵗豪航o出節(jié)點數(shù)據(jù)，，，，作三階牛頓插值多項式，計算，并估計其誤差.解首先將名為newdscg.m的程序保存為M文件，然后在MATLAB工作窗口輸入程序>>symsM,X=[-4,0,1,2];Y=[27,1,2,17];x=-2.345;[y,R,A,C,P]=newdscg(X,Y,x,M)運行后輸出插值y及其誤差限公式R，三階牛頓插值多項式P及其系數(shù)向量C，差商的矩陣A如下y=22.3211R=1323077530165133/562949953421312*M（即R=2.3503*M）A=27.00000001.0000-6.5000002.00001.00001.5000017.000015.00007.00000.9167C=0.91674.2500-4.16671.0000P=11/12*x^3+17/4*x^2-25/6*x+1&4龍貝格積分法的上機實現(xiàn)實例：clcclearsyms

x;

n=8;a=0;b=1;R=0.5*10^(-6);

T=zeros(n,n);

y=sin(x)/x;

p=subs(y,x,a);

q=subs(y,x,b);

if

a==0

p=1;

end

T(1,1)=(b-a)/2*(p+q);

for

i=2:n

f=0;

for

j=1:2^(i-2)

t=a+((2*j-1)/2^(i-1))*(b-a);

z=subs(y,x,t);

f=f+z;

end

T(i,1)=0.5*T(i-1,1)+f*(b-a)/2^(i-1);

end

for

j=2:n

n=n-1;

for

i=1:n

T(i,j)=(4^(j-1)*T(i+1,j-1)-T(i,j-1))/(4^(j-1)-1);

end

end

T=vpa(T,7)

T=eval(T);

n=8;

w=ones(1,7);

for

j=2:n

if

abs(T(1,j)-T(1,j-1))<=R

end

end

a=T(1,j)；

a=vpa(a,7)

&6

矩陣的直接LU分解法解方程組和微分方程的數(shù)值解法上機實現(xiàn)（可以用matlab自帶的結(jié)法解微分方程組）：編輯一個LU分解函數(shù)如下：function[L,U]=Lu(A)%求解線性方程組的三角分解法%A為方程組的系數(shù)矩陣%L和U為分解后的下三角和上三角矩陣[n,m]=size(A);ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!');return;end%判斷矩陣能否LU分解forii=1:nfori=1:iiforj=1:iiAA(i,j)=A(i,j);endendif(det(AA)==0)error('ThematrixcannotbedividedbyLU!')return;endend%開始計算，先賦初值L=eye(n);U=zeros(n,n);%計算U的第一行，L的第一列fori=1:nU(1,i)=A(1,i);L(i,1)=A(i,1)/U(1,1);end%計算U的第r行，L的第r列fori=2:nforj=i:nfork=1:i-1M(k)=L(i,k)*U(k,j);endU(i,j)=A(i,j)-sum(M);endforj=i+1:nfork=1:i-1M(k)=L(j,k)*U(k,i);endL(j,i)=(A(j,i)-sum(M))/U(i,i);endend然后，編輯一個通過LU分解法解線性方程組的函數(shù)如下function[L,U,x]=Lu_x(A,d)%三角分解法求解線性方程組，LU法解線性方程組Ax=LUx=d%A為方程組的系數(shù)矩陣%d為方程組的右端項%L和U為分解后的下三角和上三角矩陣%x為線性方程組的解[n,m]=size(A);ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!');return;end%判斷矩陣能否LU分解forii=1:nfori=1:iiforj=1:iiAA(i,j)=A(i,j);endendif(det(AA)==0)error('ThematrixcannotbedividedbyLU!')return;endend[L,U]=Lu(A);%直接調(diào)用自定義函數(shù)，首先將矩陣分解，A=LU%設(shè)Ly=d由于L是下三角矩陣,所以可求y(i)y(1)=d(1);fori=2:nforj=1:i-1d(i)=d(i)-L(i,j)*y(j);endy(i)=d(i);end%設(shè)Ux=y,由于U是上三角矩陣，所以可求x(i)x(n)=y(n)/U(n,n);fori=(n-1):-1:1forj=n:-1:i+1y(i)=y(i)-U(i,j)*x(j);endx(i)=y(i)/U(i,i);end然后，n=5時，調(diào)用自定義函數(shù)[L,U,x]=Lu_x(A,a)解出：x=0.9999999999896551.0000000000326090.9999999999616251.0000000000199840.999999999996114[L,U,x]=Lu_x(B,b)解出：x=0.9999999999999761.0000000000003420.9999999999988191.0000000000014770.999999999999388微分方程的數(shù)值解法上機實現(xiàn)