(完整word版)数值分析的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
for 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)*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 h>1 or h<0
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);
y=y1+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;
end
t=[0:h:T];
% 红线为解析解，'x'为数值解
plot(t,S1,'r',t,S3,'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 h>1 or h<0
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);
l2=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;
end
t=[0:h:T];
% 红线为解析解，'x'为数值解
plot(t,S1,'r',t,S3,'x')
高斯-赛德尔
function gaosisaideer
n=input('n='); %阶数
tol=input('tol='); %迭代精度
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
b(i,1)=sum(A(i,:));
end
y=zeros(n,1); %迭代解
x1=zeros(n,1); %准确解
for i=1:n
y(i,1)=0; %迭代解赋初值
x1(i,1)=1; %生成准确解
end
k=0;
while norm(y-x1)>=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);
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')
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
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 %迭代开始
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):h=') ; if h>1 or h<0
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')。