数值分析的MATLAB程序

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')