数值分析的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