數(shù)值分析的MATLAB程序

1、列主元法%選擇階數(shù)%系數(shù)矩陣A%矩陣b%解X%生成hilbert矩陣A%生成矩陣bfun ctio n lia nzhuyua n( A,b) n=input('請(qǐng)輸入 n:') A=zeros( n,n); b=zeros( n,1); X=zeros( n,1);for i=1: nfor j=1: nA(i,j)=(1/(i+j-1)； end b(i,1)=sum(A(i,:);end for i=1: n-1j=i;top=max(abs(A(i: n,j);% 列主元k=j;while abs(A(k,j)=top % 列主元所在行 k=k+1;endfor z=1

2、: n%交換主元所在行a1=A(i,z);A(i,z)=A (k, z);A(k, z)=a1;enda2=b(i,1);b(i,1)=b(k,1);b(k,1)=a2;for s=i+1:n%消去算法開(kāi)始m=A(s,j)/A(i,j);%化簡(jiǎn)為上三角矩陣A(s,j)=0;for p=i+1: nA(s,p)=A(s,p)-m*A(i,p);endb(s,1)=b(s,1)-m*b(i,1);endendX( n,1)=b( n,1)/A( n,n);% 回代開(kāi)始for i=n-1:-1:1s=0;%初始化sfor j=i+1: ns=s+A(i,j)*X(j,1);endX(i,1)=(b(

3、i,1)-s)/A(i,i);end歐拉法clcclear%歐拉法p=10;%貝塔的取值T=10;%t取值的上限y1=1;%y1的初值r1=1;%y2的初值%輸入步長(zhǎng)h的值h=i nput('歐拉法 please in put number(h=1 0.5 0.25 0.125 0.0625):h='); if h>1 or h<0breakendS1=0:T/h;S2=0:T/h;S3=0:T/h;S4=0:T/h;i=1;%迭代過(guò)程for t=0:h:TY=(exp(-t);R=(1/(p-1)*exp(-t)+(p-2)/(p-1)*exp(-p*t);y=y

4、1+h*(-y1);y1=y;r=r1+h*(y1-p*r1);r1=r;S1(i)=Y;S2(i)=R;S3(i)=y;S4(i)=r;i=i+1;endt=0:h:T;%紅線為解析解，x'為數(shù)值解plot(t,S1,'r',t,S3,'x')改進(jìn)歐拉法clcclearp=10;%貝塔的取值T=10;%t取值的上限y1=1;%y1的初值r1=1;%y2的初值%輸入步長(zhǎng)h的值h=i nput('改進(jìn)歐拉法 please in put number(h=1 0.5 0.25 0.125 0.0625):h='); if h>1 or

5、h<0breakendS1=0:T/h;S2=0:T/h;S3=0:T/h;S4=0:T/h;i=1;% 迭代過(guò)程for t=0:h:TY=(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);I2=y1+h*k1-p*(r1+h*l1);y=y1+h*(k1+k2)/2;r=r1+h*(k1+k2)/2;r1=r;y1=y;S1(i)=Y;S2(i)=R;S3(i)=y;S4(i)=r;i=i+1;endt=0:h:T;%紅線為解析解，x'為數(shù)值解plot(t,S1,&

6、#39;r',t,S3,'x')高斯-賽德?tīng)杅un ctio n gaosisaideer n=in put(' n=');tol=i nput('tol=');A=zeros( n,n); b=zeros( n,1);for i=1: nfor j=1: nA(i,j)=(1/(i+j-1)； end%階數(shù)%迭代精度%生成b向量%給Hilbert矩陣和b向量賦值b(i,1)=sum(A(i,:);endy=zeros( n,1);% 迭代解x1=zeros( n,1);% 準(zhǔn)確解for i=1: ny(i,1)=0;%迭代解賦初值x1(

7、i,1)=1;%生成準(zhǔn)確解endk=0;while norm(y-x1)>=tol %精度控制(采用自動(dòng)步數(shù)控制)k=k+1;for i=1: n%迭代開(kāi)始a1=0;a2=0;for j=1:i-1a仁 a1+A(i,j)*y(j,1);end for j=i+1: na2=a2+A(i,j)*y(j,1);endy(i,1)=(b(i,1)-a1-a2)/A(i,i);endenddisp('迭代步數(shù)k')kdisp(y)%顯示 yend最速下降法fun ctio n gaosisaideern=input('階數(shù) n=');%階數(shù)tol=input(&

8、#39;迭代精度tol=');%迭代精度eps=in put(最 速下降法 eps=');A=zeros( n,n);b=zeros(n,1);%生成 b 向量for i=1:n%給Hilbert矩陣和b向量賦值for j=1: nA(i,j)=(1/(i+j-1)；endb(i,1)=sum(A(i,:);endy=zeros( n,1);% 迭代解x1=zeros( n,1);% 準(zhǔn)確解t=zeros( n,1); r=zeros( n,1);for i=1: n y(i,1)=0; x1(i,1)=1;end%迭代解賦初值%生成準(zhǔn)確解r=b-A*y;while no rm

9、(r)>=eps;%先進(jìn)行最速下降法求得進(jìn)行賽德?tīng)柕某跏冀鈟t=(r'*r)/(r'*A*r);s1=t*r;y=y+s1;r=b-A*y;endk=0;while norm(y-x1)>=tol %精度控制(采用自動(dòng)步數(shù)控制) k=k+1;for i=1: n%迭代開(kāi)始a仁0;a2=0;for j=1:i-1a1=a1+A(i,j)*y(j,1);endfor j=i+1: na2=a2+A(i,j)*y(j,1);end y(i,1)=(b(i,1)-a1-a2)/A(i,i);endenddisp('迭代步數(shù)k') disp(k)%顯示yd

10、isp(y)四階龍格-庫(kù)塔法clcclearp=10;%貝塔的取值T=10;%t取值的上限y1=1;%y1的初值r1=1;%y2的初值%輸入步長(zhǎng)h的值h=i nput('四階龍格 please in put number(h=1 0.5 0.25 0.125 0.0625):h='); if h>1 or h<0breakendS1=0:T/h;S2=0:T/h;S3=0:T/h;S4=0:T/h;i=1;% 迭代過(guò)程for t=0:h:TY=(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);I4=y1+h*k3-p*(r1+h*l3);y=y1+h