《数值分析的MATLAB程序.docx》由会员分享，可在线阅读，更多相关《数值分析的MATLAB程序.docx（6页珍藏版）》请在三一办公上搜索。
1、数值分析的MATLAB程序列主元法 function lianzhuyuan(A,b) n=input(请输入n:) %选择阶数 A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1); %生成hilbert矩阵A end b(i,1)=sum(A(i,:); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j); %列主元 k=j; while abs(A(k,j)=top %列主元所在行 k=k+1; end f
2、or z=1:n %交换主元所在行 a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始 m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j
3、)*X(j,1); end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input(欧拉法 please input number(h=1 0.5 0.25 0.125 0.0625):h=) ; if h1 or h1 or h=tol %精度控制(采用自动步数控制) k=k+1; for i=1:n %迭代开始 a1=0; a2=0; for j=1:i-1 a1=a1+A(i,j)*y(j,1); en
4、d for j=i+1:n a2=a2+A(i,j)*y(j,1); end y(i,1)=(b(i,1)-a1-a2)/A(i,i); end end disp(迭代步数k) k disp(y) %显示y end 最速下降法 function gaosisaideer n=input(阶数n=); %阶数 tol=input(迭代精度tol=); %迭代精度 eps=input(最速下降法eps=); A=zeros(n,n); b=zeros(n,1); %生成b向量 for i=1:n %给Hilbert矩阵和b向量赋值 for j=1:n A(i,j)=(1/(i+j-1); end
5、b(i,1)=sum(A(i,:); end y=zeros(n,1); %迭代解 x1=zeros(n,1); %准确解 t=zeros(n,1); r=zeros(n,1); for i=1:n y(i,1)=0; %迭代解赋初值 x1(i,1)=1; %生成准确解 end r=b-A*y; while norm(r)=eps; %先进行最速下降法求得进行赛德尔迭代的初始解y t=(r*r)/(r*A*r); s1=t*r; y=y+s1; r=b-A*y; end k=0; while norm(y-x1)=tol %精度控制(采用自动步数控制) k=k+1; for i=1:n %迭代
6、开始 a1=0; a2=0; for j=1:i-1 a1=a1+A(i,j)*y(j,1); end for j=i+1:n a2=a2+A(i,j)*y(j,1); end y(i,1)=(b(i,1)-a1-a2)/A(i,i); end end disp(迭代步数k) disp(k) disp(y) %显示y 四阶龙格-库塔法 clc clear p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input(四阶龙格 please input number(h=1 0.5 0.25 0.125 0.0625)
7、:h=) ; if h1 or h0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t); R=(1/(p-1)*exp(-t)+(p-2)/(p-1)*exp(-p*t); k1=-y1; l1=y1-p*r1; k2=-(y1+h*k1/2); l2=y1+h*k1/2-p*(r1+h*l1/2); k3=-(y1+h*k2/2); l3=y1+h*k2/2-p*(r1+h*l2/2); k4=-(y1+h*k3); l4=y1+h*k3-p*(r1+h*l3); y=y1+h*(k1+2*k2+2*k3+k4)/6; r=r1+h*(l1+2*l2+2*l3+l4)/6; r1=r; y1=y; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=0:h:T; % 红线为解析解，x为数值解 plot(t,S1,r,t,S3,x)
