1、数值计算方法实验报告 数值计算方法实验报告叶耀 北师珠 2016.01.10&1：实验一:基本变量与运算：for 循环：s=0 for i=1:n（n为常数） s=s+i;end s函数及函数实验方法：分段函数的例子：function y = func(x)if x0 y=sin(x);else y=1-cos(x);end基本绘图命令实验：下面看如何画出这个分段函数的图形：x=-pi:0.1:pi;y=;for i=1:length(x) y(i)=func(x(i);endplot(x,y)画图的程序：x=-5:0.1:5 ;y=;for i=1:length(x) y(i)=func2(
2、x(i);endplot(x,y)grid on&2：lagrange 插值多项式实现：完善的拉格朗日插值基函数的程序：function f = language( x,y,x0 )%the Language ploynomial function % the vector x and y are knownsyms t;if (length(x)=length(y) n=length(x);else disp(x and y have diffierent dimensions !) return;end f=0;for (i=1:n) l=y(i); for(j=1:i-1) l=l*(t
3、-x(j)/(x(i)-x(j); end for(j=i+1:n) l=l*(t-x(j)/(x(i)-x(j); end f=f+l; simplify(f); if (i=n) if (nargin=3) f=subs(f,t,x0); else f=collect(f); f=vpa(f,6); end end endend&3：Newton 插值多项式的程序：function y,R,A,C,L=newdscg(X,Y,x,M)n=length(X); m=length(x);for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y;s=0.0; p=1.0;
4、q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1)/(X(i)-X(i-j+1); end q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n);for k=(n-1):-1:1C=conv(C,poly(X(k);d=length(C);C(d)=C(d)+A(k,k);end y(k)= polyval(C, z);endR=M*q1/c1;L(k,:)=poly2sym(C)&4：龙贝格积分法的上机实现：format longf=input(请输
5、入原函数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);for k=2:4sum1=0;for i=1:2(k-2)sum1=sum1+subs(sym(f),findsym(sym(f),(a+(2*i-1)*(b-a)/2(k-1);endT(k)=1/2*T(k-1)+(b-a)/(2(k-1)*sum1;endfor k=1:3S(k)=T(k+1)+1/
6、(4-1)*(T(k+1)-T(k);endfor k=1:2c(k)=S(k+1)+1/(42-1)*(S(k+1)-S(k);endR(1)=C(2)+1/(43-1)*(C(2)-C(1);k=3;while 1T(1)=T(2);T(2)=T(3);T(3)=T(4);sum2=0;for i=1:2ksum2=sum2+subs(sym(f),findsym(sym(f),(a+(2*i-1)*(b-a)/2(k+1);endT(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)
7、;C(1)=C(2);C(2)=S(3)+1/(42-1)*(S(3)-S(2);R(2)=C(2)+1/(43-1)*(C(2)-C(1);if abs(R(2)-R(1) syms M,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及其系数向量C，差商的矩阵A如下y = 22.3211R =1323077530165133/562949953421312*M（即R =2.3503*M）A= 27.0000 0 0 0 1.0000 -6.5000 0 0 2.
8、0000 1.0000 1.5000 0 17.0000 15.0000 7.0000 0.9167C = 0.9167 4.2500 -4.1667 1.0000P =11/12*x3+17/4*x2-25/6*x+1&4龙贝格积分法的上机实现实例：clcclearsymsx;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);ifa=0p=1;endT(1,1)=(b-a)/2*(p+q);fori=2:nf=0;forj=1:2(i-2)t=a+(2*j-1)/2(i-1)*(b-a);
9、z=subs(y,x,t);f=f+z;endT(i,1)=0.5*T(i-1,1)+f*(b-a)/2(i-1);endforj=2:nn=n-1;fori=1:nT(i,j)=(4(j-1)*T(i+1,j-1)-T(i,j-1)/(4(j-1)-1);endendT=vpa(T,7)T=eval(T);n=8;w=ones(1,7);forj=2:nifabs(T(1,j)-T(1,j-1)=Rendenda=T(1,j)；a=vpa(a,7)&6矩阵的直接LU分解法解方程组和微分方程的数值解法上机实现（可以用matlab自带的结法解微分方程组）：编辑一个LU分解函数如下 ： funct
10、ionL,U=Lu(A) % 求解线性方程组的三角分解法 % A为方程组的系数矩阵 %L和U为分解后的下三角和上三角矩阵 n,m=size(A); if n=m error(The rows and columns of matrix A must be equal!); return; end %判断矩阵能否LU分解 for ii=1:n for i=1:ii for j=1:ii AA(i,j)=A(i,j); end end if (det(AA)=0)error(The matrix can not be divided by LU!) return; end end %开始计算，先赋
11、初值 L=eye(n); U=zeros(n,n); %计算U的第一行，L的第一列 for i=1:n U(1,i)=A(1,i); L(i,1)=A(i,1)/U(1,1); end %计算U的第r行，L的第r列 for i=2:n for j=i:n for k=1:i-1 M(k)=L(i,k)*U(k,j); end U(i,j)=A(i,j)-sum(M); end for j=i+1:n for k=1:i-1 M(k)=L(j,k)*U(k,i); end L(j,i)=(A(j,i)-sum(M)/U(i,i); end end然后，编辑一个通过LU分解法解线性方程组的函数如下
12、 function L,U,x=Lu_x(A,d) %三角分解法求解线性方程组，LU法解线性方程组Ax=LUx=d %A为方程组的系数矩阵 %d为方程组的右端项 %L和U为分解后的下三角和上三角矩阵%x为线性方程组的解 n,m=size(A); if n=m error(The rows and columns of matrix A must be equal!); return;end %判断矩阵能否LU分解 for ii=1:n for i=1:ii for j=1:iiAA(i,j)=A(i,j); end end if (det(AA)=0) error(The matrix can
13、 not be divided by LU!) return; endend L,U=Lu(A); %直接调用自定义函数，首先将矩阵分解，A=LU %设Ly=d由于L是下三角矩阵,所以可求y(i) y(1)=d(1);for i=2:n for j=1:i-1 d(i)=d(i)-L(i,j)*y(j); end y(i)=d(i);end %设Ux=y,由于U是上三角矩阵，所以可求x(i) x(n)=y(n)/U(n,n);for i=(n-1):-1:1 for j=n:-1:i+1 y(i)=y(i)-U(i,j)*x(j); end x(i)=y(i)/U(i,i);end然后，n=5
14、时，调用自定义函数 L,U,x=Lu_x(A,a) 解出： x =0.999999999989655 1.000000000032609 0.999999999961625 1.000000000019984 0.999999999996114 L,U,x=Lu_x(B,b) 解出： x =0.999999999999976 1.000000000000342 0.999999999998819 1.000000000001477 0.999999999999388微分方程的数值解法上机实现程序：clc u=zeros(11,7); u1=zeros(11,7);pi=3.1415926; h=0.1; k=0.1; r=k/h;for i=1:11 u(i,1)=(1/8)*sin(pi*(i-1)*h) end for i=1:11 u(i,2)=(1/8)*sin(pi*(i-1)*h)+r*r*(1/8)*sin(pi*(i-1)*h)*(cos(pi*h)-1); end for j=1:7 u(1,j)=0; end for j=1:7 u(10,j)=0; end for j=2:6
