计算方法常用计算编程

X = yX1i1.向量范数:i=1IXI2、丄=max x1< i < n 'fn=yclearx=input('输入一个向量\n x=') x1=sum(abs(x))%1-范数 x2=sqrt(sum(abs(x).人2))%2-范数 xinf=max(abs(x))% 无穷-范数 p=input('多少范数?输入p=') xp=(sum(abs(x).Ap)).A(1/p)%p-范数2.矩阵范数:aij丿max(atA )ijaIIfn=max y1derta) x2=inv(D)*(L+U)*x1+inv(D)*b;t=max(abs(x1-x2)); x1=x2;endinput('回车后输出雅克比迭代法矩阵迭代的解x=') x2%%%%%%%%%%%%%%雅克比迭代法元素迭代法% t=1;while(t>derta)for i=1:ntt=0;for j=1:i-1tt=tt+a(i,j)*x1(j);endttt=0;for j=i+1:nttt=ttt+a(i,j)*x1(j);endx2(i)=(b(i)-tt-ttt)/a(i,i);endt=max(abs(x1-x2));x1=x2;endinput('回车后输出雅克比迭代法元素迭代的解%=')x2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%高斯迭代法 矩阵迭代%t=1;while(t>derta)x2=inv(D-L)*U*x1+inv(D-L)*b;t=max(abs(x1-x2));x1=x2;endinput('回车后输出高斯迭代法矩阵迭代的解x=')x2%%%%%%%%%%%%%%高斯迭代法元素迭代法%t=1;while(t>derta)for i=1:ntt=0;for j=1:i-1tt=tt+a(i,j)*x2(j);endttt=0;for j=i+1:nttt=ttt+a(i,j)*x1(j);endx2(i)=(b(i)-tt-ttt)/a(i,i);endt=max(abs(x1-x2));x1=x2;endinput('回车后输出高斯迭代法元素迭代的解%=')x2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%逐次超松弛迭代法 矩阵迭代%t=1;w=input('输入超松弛因子闪=')if w>0&w<=2w=w;else w=0.15;endwhile(t>derta)x2=inv(D-w*L)*((1-w)*D+w*U)*x1+w*inv(D-w*L)*b;t=max(abs(x1-x2));x1=x2;endinput('回车后输出超松弛迭代法 矩阵迭代的解x=')x2%%%%%%%%%%%%%%超松弛迭代法元素迭代法%t=1;while(t>derta)for i=1:ntt=0;for j=1:i-1tt=tt+a(i,j)*x2(j);endttt=0;for j=i+1:nttt=ttt+a(i,j)*x1(j);endx2(i)=w*(b(i)-tt-ttt)/a(i,i)+x1(i);endt=max(abs(x1-x2));x1=x2;endinput('回车后输出超松弛迭代法元素迭代的解%=')x2% 输入系数矩阵 a=[-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4]%输入自由项b=[1 ;1 ;1;1]%迭代初始值x1=[0 0 0 0]%输入精度要求der ta=0.000046cleara=input('任意输入一个数a=') for k=2:a-1if rem(a,k)==0sprin tf ('%d 不是素数',a) breakendendif k==a-1sprin tf('%d 是素数',a)end%%%向量范数通用程序clearx=input('输入一个向量\n x=')x1二sum(abs(x))%1-范数x2=sqr t(sum(abs(x)."2))%2-范数 xinf=max(abs(x))%无 穷-范数 p=input('多少范数?输入p=') xp=(sum(abs(x)."p))."(1/p)%p—范数 %%%矩阵范数通用程序cleara=input('输入一个矩阵\n a=') x1=max(sum(abs(a)))%1 -范数 x2=sqr t(max(eig(a'*a)))%2-范数 xinf=max(sum(abs(a)'))%无 穷-范数 xF=sqr t(sum(sum(a.“2)))%F-范数%%%高斯顺序消去法通用程序cleara=input('输入一个系数矩阵\n a=') b=input('输入一个自由项\n b=') n=length(b);for k=1:n—1for i=k+1:nl(i,k)=—a(i,k)/a(k,k);for j=k:n a(i,j)=a(i,j)+l(i,k)*a(k,j);end b(i)=b(i)+l(i,k)*b(k);endend zengguang=[a,b] x(n)=b(n)/a(n,n); for k=n-1:-1:1t=0;for j=k+1:n t=t+a(k,j)*x(j)end x(k)=(b(k)-t)/a(k,k);endx clear clcx0=input('输入初值x0=')if abs(x0)>0x0=x0;else x0=1.5 end x=sym('x');%f='x-x"3-4*x"2+10' f=input('输入迭代函数f=') t=1;pp=1; %%%%%%下面判断迭代是否收敛 g=diff(f);x=x0;gg=eval(g)if abs(gg)>1x1='迭代格式不收敛’elsederta二input('输入迭代精度derta二') if abs(derta)>0 derta=derta;else derta=0.000001; end t=1;x=x0; while(t>derta)pp=pp+1;x1=eval(f);t=abs(x1-x);x=x1;if pp>10000 % 如果不收敛,则迭代10000次后自动终止 xl二'迭代超过10000次’; break;endendendpx1%f='x-x"3-4*x"2+10'%f='1/2*(10—x飞厂(1/2)'%f='(10/x—4*x),1/2)'%f='(10/(4+x)厂(1/2)'%f='x-(x"3+4*x"2-10)/(3*x"2+8*x)'%%%%%%%%%%%%%%%%结果展示%输入初值x0=1.5 %x0 =% 1.5000% 输入迭代函数 f='x-x"3-4*x"2+10' %f =%x-x"3-4*x"2+10%输入迭代精度der ta=0.00002%derta =% 2.0000e-005%x1 =% NaN clearclcx0=input('输入初值x0=')if abs(x0)>0x0=x0;else x0=1.5endx=sym('x');%f='x-x"3-4*x"2+10'f=input('输入迭代函数f=') t=1;pp=1; %%%%%%下面判断迭代是否收敛 g=diff(f);x=x0;gg=eval(g);if abs(gg)>1x1='迭代格式不收敛’elsederta二input('输入迭代精度derta二')if abs(derta)>0 derta=derta;else derta=0.000001;endt=1;x=x0;w=1/(1-gg);% 加权因子while(t>derta)pp=pp+1;x1=w*x+(1-w)*eval(f);% 两个组合的叠加 t=abs(x1-x);x=x1;if pp>10000 % 如果不收敛,则迭代10000次后自动终止 x1=' 迭代超过10000次';break;endendendpx1 %f='x-x"3-4*x"2+10' %f='1/2*(10—x飞厂(1/2)'%f='(10/x—4*x),1/2)'%f='(10/(4+x)厂(1/2)' %f='x-(x"3+4*x"2-10)/(3*x"2+8*x)'%%%%%%%%%%%%%%%% clearclcx0=input('输入初值x0=')if abs(x0)>0x0=x0;else x0=1.5endx=sym('x');%f='x-x"3-4*x"2+10'f=input('输入迭代函数f=') t=1;pp=1; %%%%%%下面判断迭代是否收敛 g=diff(f);x=x0;gg=eval(g);if abs(gg)>1x1='迭代格式不收敛’elsederta二input('输入迭代精度derta二')if abs(derta)>0 derta=derta;else derta=0.000001;endt=1;x=x0;w=1/(1-gg);% 加权因子while(t>derta) pp=pp+1;x2=x;y=eval(f);x=y;z=eval(f);x=x2;x1=x-(y-x厂2/(z-2*y+x);%两个组合的叠加 t=abs(x1-x);x=x1;if pp>10000 % 如果不收敛,则迭代10000次后自动终止 x1=' 迭代超过10000次';break;endendendpx1%f='x-x"3-4*x"2+10'%f='l/2*(10—x飞厂(1/2)'%f='(10/x—4*x),1/2)'%f='(10/(4+x)厂(1/2)' %f='x-(x"3+4*x"2-10)/(3*x"2+8*x)' %%%%%%%%%%%%%%%%%乘幂法 clear clcA=input('输入一个n阶方阵A=') n=length(A(1,:))v0=ones(1,n)'; u0=v0/max(v0) q=1;while( q>0.00000001)v1=A*u0j=1;t=1;for i=1:nif abs(v1(i))>abs(v1(j)) j=i;endendfor i=1:nif abs(v0(i))>abs(v0(t)) t=i;endendu1=v1./v1(j) q=abs(abs(v1(j))-abs(v0(t))); u0=u1;v0=v1;endj=1;for i=1:nif abs(v1(i))>abs(v1(j)) j=i; endendlamda=v1(j)tezhengxiangliang=v1'./v1(j)%输入一个 n 阶方阵 A=[-12 3 3;3 1 -2;3 -2 7] clearclcxO=input('输入起始节点坐标xO二') h=input('输入步长h=')y=input('输入节点坐标函数值f(x)=') x2=input('输入所要计算的节点x2=') syms tn=length(y);for i=1:nx1(i)=xO+h*(i-1);endfor i=2:n f1(i,1)=(y(i)-y(i-1));end for i=2:nfor j=i+1:n f1(j,i)=(f1(j,i-1)-f1(j-1,i-1));endendf1=[y',f1]%输出带0阶差分的差分表格 %%%%%%%%%%%%%%%%%%%%%%%%%%等距节点插值函数中牛顿向前插值函数 Newton1=f1(1,1);ttt=1;for i=2:ntt=1;for j=1:i-1tt=tt*(t-j+1);endttt=ttt*(i-1);Newton1=Newton1+f1(i,i)*tt/(ttt)endfprintf('等距节点插值函数中牛顿向前插值函数为\n')expand(Newton1) t=(x2-x0)/h;p=eval(Newton1);fprintf('等距节点插值函数中牛顿向前插值函数为在所求点x2的函数值为\n')p%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%输入起始节点坐标x0=0.4%输入步长h=0.1% 输入节点坐标函数值 f(x) = [0.38942 0.47943 0.56464]%输入所要计算的节点x2=0.43251 %Newton插值函数为为:l+17/6*x-2/3*x 2+l/6*x 4T/3*x 3%New ton插值函数在所求点x2的函数值为:23 clearclcx0=input('输入起始节点坐标x0=') h=input('输入步长h=')y=input('输入节点坐标函数值f(x)=') x2=input('输入所要计算的节点x2=') syms tn=length(y);for i=1:nx1(i)=x0+h*(i-1);endfor i=2:nf1(i,1)=(y(i)-y(i-1));end for i=2:nfor j=i+1:n f1(j,i)=(f1(j,i-1)-f1(j-1,i-1));endendf1=[y',f1]% 输出带0阶差分的差分表格%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%等距节点插值函数中牛顿向后插值函数 Newton2=f1(n,1);ttt=1;for i=2:ntt=1;for j=1:i-1tt=tt*(t-j+1);endttt=ttt*(i-1);Newton2=Newton2+f1(n,i)*tt/(ttt)endfprintf('等距节点插值函数中牛顿向后插值函数为\n')expand(Newton2)t=(x2-x0)/h;p=eval(Newton2);fprintf('等距节点插值函数中牛顿向后插值函数为在所求点x2的函数值为\n')p%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%输入起始节点坐标x0=0.4%输入步长h=0.1% 输入节点坐标函数值 f(x) = [0.38942 0.47943 0.56464]%输入所要计算的节点x2=0.43251%New ton 插值函数为为:l+17/6*x-2/3*x"2+l/6*x"4T/3*x"3%New ton插值函数在所求点x2的函数值为:23 clearclcformat longsyms xf=input('输入被积函数f=')a=inpu t('输入积分区间[a b]=') %%%%%%%%%%%%%%%%%%%%%梯形积分公式 ('梯形积分公式')x=a(1); f1=eval(f);x=a(2);f2=eval(f); fun=(f1+f2)*(a(2)-a(1))/2 %%%%%%%%%%%%%%%%%%%%%辛普森积分公式('辛普森积分公式')x=a(1);f1=eval(f);x=(a(1)+a(2))/2;f2=eval(f);x=a(2);f3=eval(f); fun=(f1+4*f2+f3)*(a(2)-a(1))/6%%%%%%%%%%%%%%%%%%%%%科特斯积分公式 ('科特斯积分公式')h=(a(2)-a(1))/4; x=a(1);f1=eval(f);x=a(1)+h;f2=eval(f);x=a(1)+2*h;f3=eval(f);x=a(1)+3*h;f4=eval(f);x=a(1)+4*h;f5=eval(f);fun=(7*f1+32*f2+12*f3+32*f4+7*f5)*(a(2)-a(1))/90%实例验证P176例7.5%输入被积函数f=sqr t(x)%输入积分区间[a b] = [0.5 1] clearclca=input('输入系数矩阵a=') b=input('输入自由项b=') x1=input('输入初始值x1=') n=length(b);if sum(abs(x1))>0x1=x1;else x1=zeros(1,n);endt=2;%%%%%%%%%%%%%%%%%雅克比迭代法 while (t>0.00001)for i=1:nq=0;for j=1:i-1 q=q+a(i,j)*x1(j);endp=0;for j=i+1:n p=p+a(i,j)*x1(j);end x2(i)=(b(i)-q-p)/a(i,i);end t=max(abs(x2-x1)); x1=x2;endx2while (t>0.00001)for i=1:nq=0;for j=1:i-1 q=q+a(i,j)*x2(j);endp=0;for j=i+1:n p=p+a(i,j)*x1(j);end x2(i)=(b(i)-q-p)/a(i,i);end t=max(abs(x2-x1)); x1=x2;endx2w=input('输入松弛因子(0〈w〈2)w=‘) if w>0 w=w;else w=0.5endwhile (t>0.00001)for i=1:n q=0;for j=1:i-1q=q+a(i,j)*x2(j);end p=0;for j=i+1:np=p+a(i,j)*x1(j);end x2(i)=w*(b(i)-q-p)/a(i,i)+x1(i);endt=max(abs(x2-x1));x1=x2;endx2%%%%%%%%%%%1 1 -4]% 输入系数矩阵 a=[-4 1 1 1;1 -4 1 1 ;1 1 -4 1; 1%a =%-4111%1-411%11-41%111-4%输入自由项b=[1 1 1 1]%b =% 1 1 1 1%输入初始值x1=[0 0 0 0]%x1 =% 0 0 0 0%雅克比迭代结果%x2 =% -1.0000 -1.0000 -1.0000 -1.0000%高斯赛德尔迭代结果%x2 =% -1.0000 -1.0000 -1.0000 -1.0000%超松弛迭代结果%输入松弛因子(0〈w〈2)w=1.3%w =% 1.3000%x2 =% -1.0000 -1.0000 -1.0000 -1.0000clearclcx=input('输入节点坐标x=')y=input('输入节点坐标函数值y=')w=input('输入权重向Sw=')n=input('输入拟合最高次数n=')for i=1:n+1for j=1:n+1p(i,j) = (w.*(x.“(i—1)))*(x.”(j—1))'endendfor i=1:n+1q(i) = (w.*y)*(x."(iT))';endC=inv(p)*q'f=poly2sym(C)%%%%%%输入节点坐标x=[1 2 3 4 5]%%%输入节点坐标函数值y=[4 4.5 6 8 8.5]%%%输入权重向量w=[ 2 1 3 1 1]%%%输入拟合最高次数n=lclearclcx1=input('输入迭代初值xl二')x2=input('输入迭代初值x2=') t=1;epuc=input('输入精度要求eupc二') if epuc>0epuc=epuc;else epuc=0.000001 endwhile(t>epuc) x3=(1+x2-0.1*exp(x1))/4 x4=(xl-x「2/8)/4 t=max(abs([x3-x1,x4-x2]));x1=x3;x2=x4;endx=[x3 x4]%输入迭代初值x1=0.225%x1 =% 0.2250%输入迭代初值x2=0%x2 =%输入精度要求eupc=0.000000001clearclcx0=input('输入迭代初值xO二')x1=input('输入迭代初值x1=') t=1;epuc=input('输入精度要求eupc二') if epuc>Oepuc=epuc;else epuc=O.OOOOO1endsyms x2 x3 %x3 x4 x5fl=input('输入第一个函数(以x2 x3为变量)fl=') f2=input('输入第二个函数(以x2 x3为变量)f2=') %f3=input('输入第三个函数f3二')%f4=input('输入第四个函数f4二') f=[fl;f2]invf=jacobian(f,[x2 x3])while(t>epuc)x2=xOx3=xl xxl=[x2;x3]-inv(eval(invf))*eval(f) t=max(abs(xxl-[x2;x3]));xO=xxl(l);xl=xxl(2);endxxl%输入迭代初值x1=0%输入迭代初值x2=0%输入精度要求eupc=O.OOOOOOOO1%输入第一个函数 f1=4*x2-x3+0.1*exp(x2)T%输入第二个函数 f2=-x2+4*x3+x2"2/8 clearclcxO=input('输入迭代初值x0=')syms xf=input('输入函数f(x)=')df=diff(f)t=1;epuc=input('输入精度要求eupc二') if epuc>0epuc=epuc;else epuc=0.000001end%牛顿迭代法while(t>epuc)x=x0x1=x-eval(f)/eval(df);t=abs(x-x1);x0=x1;endx1%%%%%%%%%%%%%%%%%%%%实验数据%输入迭代初值x0=1.5%x0 =% 1.5000%输入函数 f(x)=x"3-x"2-1%f =%x"3-x"2-1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%t=1 clearclcx0=input('输入迭代初值x0=')x1=input('输入迭代初值x1=')x2=input('输入迭代初值x2=')syms xf=input('输入函数f(x)=')t=1;epuc=input('输入精度要求eupc=') if epuc>0epuc=epuc;else epuc=0.000001;end%弦割法迭代法while(t>epuc)x=x0;f0=eval(f);x=x1;f1=eval(f);x=x2;f2=eval(f);lamda3=(x2-x1)/(x1-x0);derta=1+lamda3;a=f0*lamda3“2-fl*lamda3*derta+f2*lamda3;b=f0*lamda3“2-fl*derta"2+f2*(lamda3+derta); c=f2*derta;lamda4=(-2*c)/(b+sign(b)*sqr t(b“2-4*a*c));x3=x2+lamda4*(x2-xl);t=abs(x3-x2);x0=xl;xl=x2;x2=x3;endx3%%%%%%%%%%%%%%%%%%%%实验数据%输入迭代初值x0=1.0%输入迭代初值x1=1.5%输入迭代初值x2=2.0%输入函数 f(x)=x"3+4*x"2-10%输入精度要求eupc=0.000000001%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%t=1 clearclcx0=input('输入迭代初值x0=') x1=input('输入迭代初值x1=') syms xf=input('输入函数f(x)=') t=1;epuc=input('输入精度要求eupc二') if epuc>0epuc=epuc;else epuc=0.000001;end%弦割法迭代法while(t>epuc)x=x0;f0=eval(f);x=x1;f1=eval(f);x2=x1-(f1*(x1-x0))/(f1-f0); t=abs(x2-x1);x0=x1;x1=x2;endx2%%%%%%%%%%%%%%%%%%%%实验数据%输入迭代初值x0=1.0%输入迭代初值x1=2.0%输入函数 f(x)=x"3+4*x"2-10%输入精度要求eupc=0.000000001%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%t=1 clearclcformat longsyms xf=input('输入被积函数f=')a=inpu t('输入积分区间[a b]=')n=input('输入等分份数(如果要得到辛普森积分,n应该为偶数)n=')%%%%%%%%%%%%%%%%%%%%%梯形积分公式('复化梯形积分公式')h=(a(2)-a(1))/n;t=0;x=a(1);f1=eval(f);x=a(2);f2=eval(f);for k=1:n-1x=a(1)+k*h; t=t+eval(f);endT=h*(f1+f2+2*t)/2%%%%%%%%%%%%%%%%%%%%%辛普森积分公式('辛普森积分公式') h=(a(2)-a(1))/n;x=a(1);f1=eval(f);x=a(2);f2=eval(f);t=0;for k=1:n/2x=a(1)+(2*k-1)*h; t=t+eval(f);endtt=0;for k=1:n/2-1x=a(1)+2*k*h;tt=tt+eval(f);endS=h*(f1+f2+4*t+2*tt)/6%%%%%%%%%%%%%%%%%%%%%科特斯积分公式('科特斯积分公式')h=(a(2)-a(1))/(4*n);x=a(1);f1=eval(f);x=a(2);f2=eval(f);t1=0;for k=1:n/4-1x=a(1)+(4*k+1)*h;t1=t1+eval(f);t2=0;for k=1:n/4-1 x=a(1)+(4*k+2)*h; t2=t2+eval(f);endt3=0;for k=1:n/4-1 x=a(1)+(4*k+3)*h; t3=t3+eval(f);endI=(7*f1+32*t1+12*t2+32*t3+7*f2)*h/90%实例验证P179例7.9%输入被积函数f=4/(1+x"2)%输入积分区间[a b] = [0 1] cleara=input('输入一个系数矩阵\n a=')b=input('输入一个自由项\n b=') n=length(b);for k=1:n-1for i=k+1:n l(i,k)=-a(i,k)/a(k,k);for j=k:na(i,j)=a(i,j)+l(i,k)*a(k,j); endb(i)=b(i)+l(i,k)*b(k);endendzengguang=[a,b] x(n)=b(n)/a(n,n);for k=n-1:-1:1t=0;for j=k+1:nt=t+a(k,j)*x(j)endx(k)=(b(k)-t)/a(k,k);clearclcx1=input('输入节点坐标x=')y=input('输入节点坐标函数值f(x)=') x2=input('输入所要计算的节点x2=') syms xn=length(x1);l=0;for i=1:nt1=1;for j=1:n if j==i else endendt1=t1;t1=t1*(x-x1(j));t2=1;for j=1:n if j==i else endendl=l+t1/t2*y(i);endfprintf('拉格朗日插值函数为\n')t2=t2;t2=t2*(x1(i)-x1(j));expand(l)x=x2;p=eval(l);fprintf('插值函数在所求点x2的函数值为\n')p%输入节点坐标x=[0 1 2 4]%输入节点坐标函数值f(x) = [l 9 23 3]%输入所要计算的节点x2=4%拉格朗日插值函数为:T1/4*x"3+45/4*x"2T/2*x+1%插值函数在所求点x2的函数值为clearclcformat longsyms xf=input('输入被积函数f=') a=inpu t('输入积分区间[a b]=') NN=input('输入整数(n〉2)n=‘) n=2"NN;%%%%%%%%%%%%%%%%%%%%% h=(a(2)-a(1))/2;x=a(1);f1=1;x=a(2);f2=eval(f)T(1)=(a(2)-a(1))/2*(f1+f2);N=log2(n)for k=1:Nt=0;for kk=l:2“(k-l)x=a(l) + (2*kk-l)*(a(2)-a(l))/(2*(2,k-l))); t=t+eval(f);endT(2"k)=l/2*T(2"(k-l))+h/(2"(k-l))* t;endfor k=l:N S(2"(k-l))=4*T(2"k)/3-T(2"(k-l))/3;endfor k=l:N-l C(2"(k-l))=l6*S(2"k)/l5-S(2"(k-l))/l5;endfor k=l:N-2 R(2,k-l))=64*C(22)/63-C(2,k-l))/63;endR for k=l:N+lTT(k)=T(2,k-l));for k=1:NSS(k)=S(2,k-l));endfor k=l:N-lCC(k)=C(2,k-1));endfor k=1:N-2RR(k)=R(2,k-1));end(' 梯形计算结果')TT(' 辛普森计算结果')SS(' 柯特斯计算结果')CC(' 龙贝格计算结果')RR%输入被积函数f=sin(x)/x%输入积分区间[a b] = [0 1]%输入整数n=3clearclcx1=input('输入节点坐标x=') y=input('输入节点坐标函数值f(x)=') x2=input('输入所要计算的节点x2=') syms xn=length(x1);%%%%%%%%%%%%%%%%%%%%差商的求法for i=2:n f1(i,1)=(y(i)-y(i-1))/(x1(i)-x1(i-1));endfor i=2:nfor j=i+1:n f1(j,i)=(f1(j,i-1)-f1(j-1,i-1))/(x1(j)-x1(j-i));endendf1=[y',f1]% 输出带0阶差商的差商表格%%%%%%%%%%%%%%%%%%%%%%%%%%New ton 插值函数Newton=f1(1,1);for i=2:ntt=1;for j=1:i-1 tt=tt*(x-x1(j));end Newton=Newton+f1(i,i)*tt;endfprin tf('New ton 插值函数为\n') expand(Newton)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%x=x2;p=eval(Newton);fprin tf('New ton插值函数在所求点x2的函数值为\n')p%输入节点坐标x=[-l 0 1 2 3]%输入节点坐标函数值f(x) = [-2 1 3 4 8]%输入所要计算的节点x2=4%New ton 插值函数为为:l+17/6*x-2/3*x"2+l/6*x"4T/3*x"3%New ton插值函数在所求点x2的函数值为:23 clearclca=input('输入对称正定矩阵a=')b=input('输入自由项b=') n=length(a(:,1));for k=1:nif (det(a(1:k,1:k))<=0)input('矩阵不是正定矩阵,请重新运行程序') endend%%%%%%%%%%%分 解人丸札’for i=1:nt=0;for s=1:i-1t二t+L(i,s厂2;endL(i,i)=sqrt(a(i,i)-t);for k=i+1:ntt=0;for s=1:i-1tt=tt+L(i,s)*L(k,s); endL(k,i)=(a(k,i)-tt)/L(i,i);endend%%%%%%%%%%%%%%%分解AX=b为Ly=b Lx=y%%%%%%%%%%求yfor i=1:nttt=0;for k=1:i-1 ttt=ttt+L(i,k)*y(k);endy(i)=(b(i)-ttt)/L(i,i);endfor i=n:-1:1tttt=0;for k=i+1:ntttt=tttt+L(k,i)*x(k);endx(i)=(y(i)-tttt)/L(i,i);endx clearclcb=input('输入严格三对角占优矩阵主对角线的元>b=')a=input('输入严格三对角占优矩阵主对角线下面一斜行的元素&=') c=input('输入严格三对角占优矩阵主对角线上面一斜行的元素c=') d=input('输入自由项d=')n=length(d);u(1)=b(1);for i=2:nL(i)=a(i-1)/u(i-1);u(i)=b(i)-L(i)*c(i-1);u y(1)=d(1); for i=2:n y(i)=d(i)-L(i)*y(i-1);endyx(n)=y(n)/u(n);for i=n-1:-1:1x(i)=(y(i)-c(i)*x(i+1))/u(i);endx%%%%运行结果显示%%输入严格三对角占优矩阵主对角线的元>b=[3 4 2 1]%%b =%%3 4 2 1%%输入严格三对角占优矩阵主对角线下面一斜行的元>a=[1 1 -1]%%a =%% 1 1 -1%%输入严格三对角占优矩阵主对角线上面一斜行的元>c=[2 2 1]%%c =%% 2 2 1%%输入自由项b=[5 7 2 -2]%%d =%% 5 7 2 -2%% L =%% 0 0.3333 0.3000 -0.7143%% u =%%3.00003.33331.40001.7143%% y =%%5.00005.33330.4000-1.7143%% x =%%1.00001.00001.0000-1.0000%%%% cleara=[1,2,3;1,3,5;1,3,6] n=length(a(1,:));for j=1:nU(1,j)=a(1,j);L(j,1)=a(j,1)/U(1,1);endfor k=1:nfor j=1:nt=0;for m=1:k-1t=t+L(k,m)*U(m,j);endU(k,j)=a(k,j)-t;endfor i=k+1:nt=0;for m=1:k-1t=t+L(i,m)*U(m,k);L(i,k)=a(i,k)-t;endendfor i=1:nL(i,i)=1;endLU%雅克比求特征值clearclca=input('输入一个n阶对称方阵A=')n=length(a(1,:))q=1;Q=eye(n);while (q>0.00001)i=1;j=2;for t=1:n-1for s=t+1:nif abs(a(t,s))>abs(a(i,j))i=t;j=s;endendendija(i,j)if a(i,i)==a(j,j) cosk=cos(pi/4);sink=sign(a(i,i))*cosk; elsed=(a(i,i)-a(j,j))/(-2*a(i,j)); t=sign(d)/(abs(d)+sqr t(d"2+l));cosk=l/sqrt(t"2+1); sink=t*cosk;endcosksinkp=eye(n);p(i,i)=cosk;p(j,j)=cosk;p(i,j)=sink;p(j,i)=-sink;p;a1=p'*a*pQ=Q*pq=abs(sum(sum(a. 2))-sum(diag(a. 2)));a=a1;enda1tezhengzhi=diag(a1)Q%%%%%实例验证书本P165例6.6结果%输入一个 n 阶对称方阵 A=[3.5 -6 5;-6 8.5 -9;5 -9 8.5]%书上例题有误,验证可以运用[p,v]=eig([3.5 -6 5;-6 8.5 -9;5 -9 8.5]察看并比对结果。