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')
相关推荐: