数值代数课设(99分)--基于Jacobi迭代，GS迭代，SOR迭代对泊松方程的求解[matlab](下)
2020-07-19 04:54
811 查看
附源代码–matlab
1.U1.m
%计算一维泊松方程的解的三种*2的方法
clear,clc;
n=input('请输入一维泊松方程的维数: '); %m为矩阵的维数（一维泊松方程的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
A=2*eye(n)+diag(a,1)+diag(a,-1);
b=zeros(n,1);
a1=input('x=0时的初值条件，U(0)=: '); %x=0时的初值条件
a2=input('x=1时的初值条件，U(1)=: '); %x=1时的初值条件
syms x;
f(x)=x^2+x+1+exp(x); %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
b(1)=f(X(1))*h^2+a1; %算出AX=b中的b
for i=2:n-1
b(i)=f(X(i))*h^2;
end
b(n)=f(X(n))*h^2+a2;
need_1 = input('请问使用什么迭代方式: '); %选择1.Jacobi iteration，2.the G-S iteration, 3.the SOR iteration
need_2 = input('请问使用该迭代的那种计算方法: '); %1.直接使用迭代矩阵G。2.用矩阵元素写出来的公式
if need_1==1
if need_2==1
[t,k,Uj,eps]=Jacobi(A,b);
plot(eps)
xlabel('迭代步数');
ylabel('误差的范数');
title('Jacobi1');
grid on
else
[t,k,Uj,eps]=Jacobi1(A,b);
plot(eps)
xlabel('迭代步数');
ylabel('误差的范数');
title('Jacobi2');
grid on
end
elseif need_1==2
if need_2==1
[t,k,Uj,eps]=GS(A,b);
plot(eps)
xlabel('迭代步数');
ylabel('误差的范数');
title('GS1');
grid on
else
[t,k,Uj,eps]=GS1(A,b);
plot(eps)
xlabel('迭代步数');
ylabel('误差的范数');
title('GS2');
grid on
end
elseif need_1==3
w = input('请输入松弛变量w: '); %松弛变量w
if need_2==1
[t,k,Uj,eps]=SOR(A,b,w);
plot(eps)
xlabel('迭代步数');
ylabel('误差的范数');
title('SOR1');
grid on
else
[t,k,Uj,eps]=SOR1(A,b,w);
plot(eps)
xlabel('迭代步数');
ylabel('误差的范数');
title('SOR2');
grid on
end
end
fprintf('迭代步数为:%d\n',k);
fprintf('cpu时间为:%f\n',t);
tt = linspace(0,1,n+2);
UU=[a1,Uj',a2];
aa1 = polyfit(tt,UU,2);
aa2 = polyfit(tt,UU,3);
x1=linspace(0,1,n^2);
y1=aa1(1)*x1.^2+aa1(2)*x1+aa1(3);
y2=aa2(1)*x1.^3+aa2(2)*x1.^2+aa2(3)*x1+aa2(4);
figure;
plot(tt,UU,'k*',x1,y1,'r-',x1,y2,'b-');
legend('原始数据','二次拟合','三次拟合')
fprintf('泊松方程的解为: ');
disp(poly2str(aa2,'x'));
2.Jacobi.m
function [t,k,Uj,eps]=Jacobi(A,b)
%Jacobi迭代
%x(k+1) = D^(-1)*(L+U)*x(k)+D^(-1)*b,k=1,2,…
%G=D^(-1)*(L+U)
tic;
k = 1; %记录迭代步数
D=diag(diag(A));
L=-tril(A)+D;
U=-triu(A)+D;
M=D;
N=L+U;
G=M\N;
a=zeros(length(b),1);
U(:,1)=G*a+M\b; %取x(0)=zeros(length(b),1);
epss=0.00001; %控制精度
eps=zeros(length(b),1);
eps(k)=norm(A*U(:,1)-b);
while (eps(k) >=epss)&&(k<1000)
U(:,k+1)=G*U(:,k)+M\b;
k = k + 1;
eps(k)=norm(A*U(:,k)-b);
end
t = toc;
Uj=U(:,k);
3.Jacobi1.m
function [t,k,Uj,eps]=Jacobi1(A,b)
%Jacobi迭代
%x(k+1) = D^(-1)*(L+U)*x(k)+D^(-1)*b,k=1,2,…
%G=D^(-1)*(L+U)
tic;
epss = 0.00001;
X(:,1) = (linspace(0,0,length(b)))';
k = 1;
count = 1;
eps(count) = norm(A*X(:,1)-b);
while eps(count) >= epss
for i = 1:length(b)
s = 0;
for j = 1:length(b)
s = s+A(i,j)*X(j,k);
end
s = s - A(i,i)*X(i,k);
X(i,k+1) = 1/A(i,i)*(b(i)-s);
end
k = k+1;
count = count + 1;
eps(count) = norm(A*X(:,k)-b);
end
t = toc;
Uj=X(:,k);
4.GS.m
function [t,k,Uj,eps]=GS(A,b)
%GS迭代
%x(k+1) = (D-L)^(-1)*U*x(k)+(D-L)^(-1)*b,k=1,2,…
%G=(D-L)^(-1)*U
tic;
k = 1; %记录迭代步数
D=diag(diag(A));
L=-tril(A)+D;
U=-triu(A)+D;
M=D-L;
N=U;
G=M\N;
a=zeros(length(b),1);
U(:,1)=G*a+M\b; %取x(0)=zeros(length(b),1);
epss=0.00001; %控制精度
eps(k)=norm(A*U(:,1)-b);
while (eps(k) >=epss)&&(k<1000)
U(:,k+1)=G*U(:,k)+M\b;
k = k + 1;
eps(k)=norm(A*U(:,k)-b);
end
t = toc;
Uj=U(:,k);
5.GS1.m
function [t,k,Uj,eps]=GS1(A,b)
%GS迭代
%x(k+1) = (D-L)^(-1)*U*x(k)+(D-L)^(-1)*b,k=1,2,…
%G=(D-L)^(-1)*U
tic;
epss = 0.00001;
X(:,1) = (linspace(0,0,length(b)))';
k = 1;
count = 1;
eps(count) = norm(A*X(:,1)-b);
while eps(count) >= epss
for i = 1:length(b)
s1 = 0;
s2 = 0;
for j = 1:length(b)
if j <=i-1
s1 = s1 + A(i,j)*X(j,k+1);
else
s2 = s2 + A(i,j)*X(j,k);
end
end
s2 = s2 - A(i,i)*X(i,k);
s = s1 + s2;
X(i,k+1) = 1/A(i,i)*(b(i)-s);
end
k = k+1;
count = count + 1;
eps(count) = norm(A*X(:,k)-b);
end
t = toc;
Uj=X(:,k);
6.SOR.m
function [t,k,Uj,eps]=SOR(A,b,w)
%GOR迭代
%x(k+1) = (1/w*D-L)^(-1)*((1-w)/w*D+U)*x(k)+(1/w*D-L)^(-1)*b,k=1,2,…
%G=(1/w*D-L)^(-1)*((1-w)/w*D+U)
tic;
k = 1; %记录迭代步数
D=diag(diag(A));
L=-tril(A)+D;
U=-triu(A)+D;
M=1/w*D-L;
N=(1-w)/w*D+U;
G=M\N;
a=zeros(length(b),1);
U(:,1)=G*a+M\b; %取x(0)=zeros(length(b),1);
epss=0.00001; %控制精度
eps(k)=norm(A*U(:,1)-b);
while (eps(k) >=epss)&&(k<1000)
U(:,k+1)=G*U(:,k)+M\b;
k = k + 1;
eps(k)=norm(A*U(:,k)-b);
end
t = toc;
Uj=U(:,k);
7.SOR1.m
function [t,k,Uj,eps]=SOR1(A,b,w)
%GOR迭代
%x(k+1) = (1/w*D-L)^(-1)*((1-w)/w*D+U)*x(k)+(1/w*D-L)^(-1)*b,k=1,2,…
%G=(1/w*D-L)^(-1)*((1-w)/w*D+U)
tic;
epss = 0.00001;
X(:,1) = (linspace(0,0,length(b)))';
k = 1;
count = 1;
eps(count) = norm(A*X(:,1)-b);
while eps(count) >= epss
for i = 1:length(b)
s1 = 0;
s2 = 0;
for j = 1:length(b)
if j <=i-1
s1 = s1 + A(i,j)*X(j,k+1);
else
s2 = s2 + A(i,j)*X(j,k);
end
end
s2 = s2 - A(i,i)*X(i,k);
s = s1 + s2;
X(i,k+1) = (1-w)*X(i,k)+w/A(i,i)*(b(i)-s);
end
k = k+1;
count = count + 1;
eps(count) = norm(A*X(:,k)-b);
end
t = toc;
Uj=X(:,k);
8.Test1.m
%一维泊松方程，三种迭代第一种方式实现对比图
clear,clc;
n=input('请输入一维泊松方程的维数: '); %m为矩阵的维数（一维泊松方程的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1);
A=2*eye(n)+diag(a,1)+diag(a,-1); %构造一维泊松方程的系数矩阵，记为A
b=zeros(n,1);
a1=input('x=0时的初值条件，U(0)=: '); %x=0时的初值条件
a2=input('x=1时的初值条件，U(1)=: '); %x=1时的初值条件
w = input('请输入松弛变量w: '); %松弛变量w
syms x;
f(x)=x^2+x+1+exp(x); %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
b(1)=f(X(1))*h^2+a1; %算出AX=b中的b
for i=2:n-1
b(i)=f(X(i))*h^2;
end
b(n)=f(X(n))*h^2+a2;
[tj,kj,Ujj,epsj]=Jacobi(A,b); %Jacobi:tj是cpu时间,kj是迭代次数,Ujj是A*X=b的解,epsj是误差的范数
t1=1:kj;
[tg,kg,Ujg,epsg]=GS(A,b); %GS:tg是cpu时间,kg是迭代次数,Ujg是A*X=b的解,epsg是误差的范数
t2=1:kg;
[ts,ks,Ujs,epss]=SOR(A,b,w); %SOR:ts是cpu时间,ks是迭代次数,Ujj是A*X=b的解,epss是误差的范数
t3=1:ks;
plot(t1,epsj,'r*-',t2,epsg,'b*-',t3,epss,'y*-')
grid on
xlabel('迭代步数');
ylabel('误差的范数');
title('三种迭代第一种方式实现对比图');
legend('Jacobi','GS','SOR')
fprintf('Jacobi迭代步数为:%d\n',kj);
fprintf('Jacobi迭代cpu时间为:%f\n',tj);
fprintf('GS迭代步数为:%d\n',kg);
fprintf('GS迭代cpu时间为:%f\n',tg);
fprintf('SOR迭代步数为:%d\n',ks);
fprintf('SOR迭代cpu时间为:%f\n',ts);
9.Test2.m
%一维泊松方程，三种迭代第二种方式实现对比图
clear,clc;
n=input('请输入一维泊松方程的维数: '); %m为矩阵的维数（一维泊松方程的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
A=2*eye(n)+diag(a,1)+diag(a,-1);
b=zeros(n,1);
a1=input('x=0时的初值条件，U(0)=: '); %x=0时的初值条件
a2=input('x=1时的初值条件，U(1)=: '); %x=1时的初值条件
w = input('请输入松弛变量w: '); %松弛变量w
syms x;
f(x)=x^2+x+1+exp(x); %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
b(1)=f(X(1))*h^2+a1; %算出AX=b中的b
for i=2:n-1
b(i)=f(X(i))*h^2;
end
b(n)=f(X(n))*h^2+a2;
[tj,kj,Ujj,epsj]=Jacobi1(A,b); %Jacobi1:tj是cpu时间,kj是迭代次数,Ujj是A*X=b的解,epsj是误差的范数
t1=1:kj;
[tg,kg,Ujg,epsg]=GS1(A,b); %GS1:tg是cpu时间,kg是迭代次数,Ujg是A*X=b的解,epsg是误差的范数
t2=1:kg;
[ts,ks,Ujs,epss]=SOR1(A,b,w); %SOR1:ts是cpu时间,ks是迭代次数,Ujj是A*X=b的解,epss是误差的范数
t3=1:ks;
plot(t1,epsj,'r*-',t2,epsg,'b*-',t3,epss,'y*-')
grid on
xlabel('迭代步数');
ylabel('误差的范数');
title('三种迭代第二种方式实现对比图');
legend('Jacobi','GS','SOR')
fprintf('Jacobi迭代步数为:%d\n',kj);
fprintf('Jacobi迭代cpu时间为:%f\n',tj);
fprintf('GS迭代步数为:%d\n',kg);
fprintf('GS迭代cpu时间为:%f\n',tg);
fprintf('SOR迭代步数为:%d\n',ks);
fprintf('SOR迭代cpu时间为:%f\n',ts);
10.Test3.m
%一维泊松方程，Jacobi迭代实现的两种方式对比
clear,clc;
n=input('请输入一维泊松方程的维数: '); %m为矩阵的维数（一维泊松方程的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
A=2*eye(n)+diag(a,1)+diag(a,-1);
b=zeros(n,1);
a1=input('x=0时的初值条件，U(0)=: '); %x=0时的初值条件
a2=input('x=1时的初值条件，U(1)=: '); %x=1时的初值条件
syms x;
f(x)=x^2+x+1+exp(x); %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
b(1)=f(X(1))*h^2+a1; %算出AX=b中的b
for i=2:n-1
b(i)=f(X(i))*h^2;
end
b(n)=f(X(n))*h^2+a2;
[tj1,kj1,Ujj1,epsj1]=Jacobi(A,b); %Jacobi:tj是cpu时间,kj是迭代次数,Ujj是A*X=b的解,epsj是误差的范数
t1=1:kj1;
[tj2,kj2,Ujj2,epsj2]=Jacobi1(A,b); %Jacobi1:tj是cpu时间,kj是迭代次数,Ujj是A*X=b的解,epsj是误差的范数
t2=1:kj2;
plot(t1,epsj1,'r*-',t2,epsj2,'b*-')
grid on
xlabel('迭代步数');
ylabel('误差的范数');
title('Jacobi迭代对比图');
legend('Jacobi','Jacobi1')
fprintf('第一种Jacobi迭代步数为:%d\n',kj1);
fprintf('第一种Jacobi迭代cpu时间为:%f\n',tj1);
fprintf('第二种Jacobi迭代步数为:%d\n',kj2);
fprintf('第二种Jacobi迭代cpu时间为:%f\n',tj2);
11.Test4.m
%一维泊松方程，GS迭代实现的两种方式对比
clear,clc;
n=input('请输入一维泊松方程的维数: '); %m为矩阵的维数（一维泊松方程的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
A=2*eye(n)+diag(a,1)+diag(a,-1);
b=zeros(n,1);
a1=input('x=0时的初值条件，U(0)=: '); %x=0时的初值条件
a2=input('x=1时的初值条件，U(1)=: '); %x=1时的初值条件
syms x;
f(x)=x^2+x+1+exp(x); %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
b(1)=f(X(1))*h^2+a1; %算出AX=b中的b
for i=2:n-1
b(i)=f(X(i))*h^2;
end
b(n)=f(X(n))*h^2+a2;
[tg1,kg1,Ujg1,epsg1]=GS(A,b); %GS:tg是cpu时间,kg是迭代次数,Ujg是A*X=b的解,epsg是误差的范数
t1=1:kg1;
[tg2,kg2,Ujg2,epsg2]=GS1(A,b); %GS1:tg是cpu时间,kg是迭代次数,Ujg是A*X=b的解,epsg是误差的范数
t2=1:kg2;
plot(t1,epsg1,'r*-',t2,epsg2,'b*-')
grid on
xlabel('迭代步数');
ylabel('误差的范数');
title('GS迭代对比图');
legend('GS','GS1')
fprintf('第一种GS迭代步数为:%d\n',kg1);
fprintf('第一种GS迭代cpu时间为:%f\n',tg1);
fprintf('第二种GS迭代步数为:%d\n',kg2);
fprintf('第二种GS迭代cpu时间为:%f\n',tg2);
12.Test5.m
%一维泊松方程，SOR迭代实现的两种方式对比
clear,clc;
n=input('请输入一维泊松方程的维数: '); %m为矩阵的维数（一维泊松方程的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
A=2*eye(n)+diag(a,1)+diag(a,-1);
b=zeros(n,1);
a1=input('x=0时的初值条件，U(0)=: '); %x=0时的初值条件
a2=input('x=1时的初值条件，U(1)=: '); %x=1时的初值条件
w = input('请输入松弛变量w: '); %松弛变量w
syms x;
f(x)=x^2+x+1+exp(x); %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
b(1)=f(X(1))*h^2+a1; %算出AX=b中的b
for i=2:n-1
b(i)=f(X(i))*h^2;
end
b(n)=f(X(n))*h^2+a2;
[ts1,ks1,Ujs1,epss1]=SOR(A,b,w); %SOR:ts是cpu时间,ks是迭代次数,Ujj是A*X=b的解,epss是误差的范数
t1=1:ks1;
[ts2,ks2,Ujs2,epss2]=SOR1(A,b,w); %SOR1:ts是cpu时间,ks是迭代次数,Ujj是A*X=b的解,epss是误差的范数
t2=1:ks2;
plot(t1,epss1,'r*-',t2,epss2,'b*-')
grid on
xlabel('迭代步数');
ylabel('误差的范数');
title('SOR迭代对比图');
legend('SOR','SOR1')
fprintf('第一种SOR迭代步数为:%d\n',ks1);
fprintf('第一种SOR迭代cpu时间为:%f\n',ts1);
fprintf('第二种SOR迭代步数为:%d\n',ks2);
fprintf('第二种SOR迭代cpu时间为:%f\n',ts2);
13.Test6.m
%一维泊松方程，SOR迭代实现不同w对比
clear,clc;
n=input('请输入一维泊松方程的维数: '); %m为矩阵的维数（一维泊松方程的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
A=2*eye(n)+diag(a,1)+diag(a,-1);
b=zeros(n,1);
a1=input('x=0时的初值条件，U(0)=: '); %x=0时的初值条件
a2=input('x=1时的初值条件，U(1)=: '); %x=1时的初值条件
syms x;
f(x)=x^2+x+1+exp(x); %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
b(1)=f(X(1))*h^2+a1; %算出AX=b中的b
for i=2:n-1
b(i)=f(X(i))*h^2;
end
b(n)=f(X(n))*h^2+a2;
ts1=zeros(1,19);
ks1=zeros(1,19);
Ujs1=zeros(n,19);
k=1;
for w = 0.1:0.1:1.9
[ts1(k),ks1(k),Ujs1(:,k),epss1(1:ks1(k),k)]=SOR(A,b,w); %SOR:ts是cpu时间,ks是迭代次数,Ujj是A*X=b的解,epss是误差的范数
t1=1:ks1(k);
epss=epss1(1:ks1(k),k);
plot(t1,epss)
xlabel('迭代步数');
ylabel('误差的范数');
title('SOR迭代');
fprintf('w=%f时SOR迭代步数为:%d\n',w,ks1(k));
fprintf('w=%f时SOR迭代cpu时间为:%f\n',w,ts1(k));
k=k+1;
figure
end
w=0.1:0.1:1.9;
a=polyfit(w,ks1,5);
fprintf('w和迭代次数的关系: ');
p=poly2str(a,'w');
disp(p);
x1=0.1:0.02:1.9;
y1=a(1)*x1.^5+a(2)*x1.^4+a(3)*x1.^3+a(4)*x1.^2+a(5)*x1+a(6);
plot(w,ks1,'k*',x1,y1,'r-');
legend('原始数据','五次拟合')
14.U2.m
%计算二维泊松方程的解的三种*2的方法
clear,clc;
n=input('请输入二维泊松方程对x,y的分割细密程度（矩阵的维数的平方根）: '); %m为矩阵的维数（二维泊松方程对应的一维的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
T=2*eye(n)+diag(a,1)+diag(a,-1);
I=eye(n);
A=kron(I,T)+kron(T,I);
B=zeros(n,n);
syms x y;
a1(x,y)=x^2+y^2+x; %初值条件
f(x,y)=x^2+x+1+exp(x)+y; %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
Y=linspace(1/(n+1),1-1/(n+1),n); %分割
for i=1:n %构造b
B(i,1)=f(X(i),Y(i))*h^2+a1(X(i),0);
B(1,i)=f(X(i),Y(i))*h^2+a1(0,Y(i));
B(i,n)=f(X(i),Y(i))*h^2+a1(X(i),1);
B(1,n)=f(X(i),Y(i))*h^2+a1(1,Y(i));
end
for i = 2:n-1
for j = 2:n-1
B(i,j)=f(X(i),Y(i))*h^2;
end
end
b=reshape(B,n^2,1);
need_1 = input('请问使用什么迭代方式: '); %选择Jacobi iteration，the G-S iteration, and the SOR iteration
need_2 = input('请问使用该迭代的那种计算方法: '); %1.直接使用迭代矩阵G。2.用矩阵元素写出来的公式
if need_1==1
if need_2==1
[t,k,Uj,eps]=Jacobi(A,b);
plot(eps)
xlabel('迭代步数');
ylabel('误差的范数');
title('Jacobi1');
grid on
else
[t,k,Uj,eps]=Jacobi1(A,b);
plot(eps)
xlabel('迭代步数');
ylabel('误差的范数');
title('Jacobi2');
grid on
end
elseif need_1==2
if need_2==1
[t,k,Uj,eps]=GS(A,b);
plot(eps)
xlabel('迭代步数');
ylabel('误差的范数');
title('GS1');
grid on
else
[t,k,Uj,eps]=GS1(A,b);
plot(eps)
xlabel('迭代步数');
ylabel('误差的范数');
title('GS2');
grid on
end
elseif need_1==3
w = input('请输入松弛变量w: '); %松弛变量w
if need_2==1
[t,k,Uj,eps]=SOR(A,b,w);
plot(eps)
xlabel('迭代步数');
ylabel('误差的范数');
title('SOR1');
grid on
else
[t,k,Uj,eps]=SOR1(A,b,w);
plot(eps)
xlabel('迭代步数');
ylabel('误差的范数');
title('SOR2');
grid on
end
end
fprintf('迭代步数为:%d\n',k);
fprintf('cpu时间为:%f\n',t);
fprintf('泊松方程的解为: ');
Uj
15.Test11.m
%计算二维泊松方程的解的三种*2的方法
clear,clc;
n=input('请输入二维泊松方程对x,y的分割细密程度（矩阵的维数的平方根）: '); %m为矩阵的维数（二维泊松方程对应的一维的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
T=2*eye(n)+diag(a,1)+diag(a,-1);
I=eye(n);
A=kron(I,T)+kron(T,I);
w = input('请输入松弛变量w: '); %松弛变量w
B=zeros(n,n);
syms x y;
a1(x,y)=x^2+y^2+x; %初值条件
f(x,y)=x^2+x+1+exp(x)+y; %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
Y=linspace(1/(n+1),1-1/(n+1),n); %分割
for i=1:n %构造b
B(i,1)=f(X(i),Y(i))*h^2+a1(X(i),0);
B(1,i)=f(X(i),Y(i))*h^2+a1(0,Y(i));
B(i,n)=f(X(i),Y(i))*h^2+a1(X(i),1);
B(1,n)=f(X(i),Y(i))*h^2+a1(1,Y(i));
end
for i = 2:n-1
for j = 2:n-1
B(i,j)=f(X(i),Y(i))*h^2;
end
end
b=reshape(B,n^2,1);
[tj,kj,Ujj,epsj]=Jacobi(A,b); %Jacobi:tj是cpu时间,kj是迭代次数,Ujj是A*X=b的解,epsj是误差的范数
t1=1:kj;
[tg,kg,Ujg,epsg]=GS(A,b); %GS:tg是cpu时间,kg是迭代次数,Ujg是A*X=b的解,epsg是误差的范数
t2=1:kg;
[ts,ks,Ujs,epss]=SOR(A,b,w); %SOR:ts是cpu时间,ks是迭代次数,Ujj是A*X=b的解,epss是误差的范数
t3=1:ks;
plot(t1,epsj,'r*-',t2,epsg,'b*-',t3,epss,'y*-')
grid on
xlabel('迭代步数');
ylabel('误差的范数');
title('三种迭代第一种方式实现对比图');
legend('Jacobi','GS','SOR')
fprintf('Jacobi迭代步数为:%d\n',kj);
fprintf('Jacobi迭代cpu时间为:%f\n',tj);
fprintf('GS迭代步数为:%d\n',kg);
fprintf('GS迭代cpu时间为:%f\n',tg);
fprintf('SOR迭代步数为:%d\n',ks);
fprintf('SOR迭代cpu时间为:%f\n',ts);
16.Test22.m
%二维泊松方程，三种迭代第二种方式实现对比图
clear,clc;
n=input('请输入二维泊松方程对x,y的分割细密程度（矩阵的维数的平方根）: '); %m为矩阵的维数（二维泊松方程对应的一维的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
T=2*eye(n)+diag(a,1)+diag(a,-1);
I=eye(n);
A=kron(I,T)+kron(T,I);
w = input('请输入松弛变量w: '); %松弛变量w
B=zeros(n,n);
syms x y;
a1(x,y)=x^2+y^2+x; %初值条件
f(x,y)=x^2+x+1+exp(x)+y; %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
Y=linspace(1/(n+1),1-1/(n+1),n); %分割
for i=1:n %构造b
B(i,1)=f(X(i),Y(i))*h^2+a1(X(i),0);
B(1,i)=f(X(i),Y(i))*h^2+a1(0,Y(i));
B(i,n)=f(X(i),Y(i))*h^2+a1(X(i),1);
B(1,n)=f(X(i),Y(i))*h^2+a1(1,Y(i));
end
for i = 2:n-1
for j = 2:n-1
B(i,j)=f(X(i),Y(i))*h^2;
end
end
b=reshape(B,n^2,1);
[tj,kj,Ujj,epsj]=Jacobi1(A,b); %Jacobi1:tj是cpu时间,kj是迭代次数,Ujj是A*X=b的解,epsj是误差的范数
t1=1:kj;
[tg,kg,Ujg,epsg]=GS1(A,b); %GS1:tg是cpu时间,kg是迭代次数,Ujg是A*X=b的解,epsg是误差的范数
t2=1:kg;
[ts,ks,Ujs,epss]=SOR1(A,b,w); %SOR1:ts是cpu时间,ks是迭代次数,Ujj是A*X=b的解,epss是误差的范数
t3=1:ks;
plot(t1,epsj,'r*-',t2,epsg,'b*-',t3,epss,'y*-')
grid on
xlabel('迭代步数');
ylabel('误差的范数');
title('三种迭代第二种方式实现对比图');
legend('Jacobi','GS','SOR')
fprintf('Jacobi迭代步数为:%d\n',kj);
fprintf('Jacobi迭代cpu时间为:%f\n',tj);
fprintf('GS迭代步数为:%d\n',kg);
fprintf('GS迭代cpu时间为:%f\n',tg);
fprintf('SOR迭代步数为:%d\n',ks);
fprintf('SOR迭代cpu时间为:%f\n',ts);
17.Test33.m
%二维泊松方程，Jacobi迭代实现的两种方式对比
clear,clc;
n=input('请输入二维泊松方程对x,y的分割细密程度（矩阵的维数的平方根）: '); %m为矩阵的维数（二维泊松方程对应的一维的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
T=2*eye(n)+diag(a,1)+diag(a,-1);
I=eye(n);
A=kron(I,T)+kron(T,I);
B=zeros(n,n);
syms x y;
a1(x,y)=x^2+y^2+x; %初值条件
f(x,y)=x^2+x+1+exp(x)+y; %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
Y=linspace(1/(n+1),1-1/(n+1),n); %分割
for i=1:n %构造b
B(i,1)=f(X(i),Y(i))*h^2+a1(X(i),0);
B(1,i)=f(X(i),Y(i))*h^2+a1(0,Y(i));
B(i,n)=f(X(i),Y(i))*h^2+a1(X(i),1);
B(1,n)=f(X(i),Y(i))*h^2+a1(1,Y(i));
end
for i = 2:n-1
for j = 2:n-1
B(i,j)=f(X(i),Y(i))*h^2;
end
end
b=reshape(B,n^2,1);
[tj1,kj1,Ujj1,epsj1]=Jacobi(A,b); %Jacobi:tj是cpu时间,kj是迭代次数,Ujj是A*X=b的解,epsj是误差的范数
t1=1:kj1;
[tj2,kj2,Ujj2,epsj2]=Jacobi1(A,b); %Jacobi1:tj是cpu时间,kj是迭代次数,Ujj是A*X=b的解,epsj是误差的范数
t2=1:kj2;
plot(t1,epsj1,'r*-',t2,epsj2,'b*-')
grid o
1c8e8
n
xlabel('迭代步数');
ylabel('误差的范数');
title('Jacobi迭代对比图');
legend('Jacobi','Jacobi1')
fprintf('第一种Jacobi迭代步数为:%d\n',kj1);
fprintf('第一种Jacobi迭代cpu时间为:%f\n',tj1);
fprintf('第二种Jacobi迭代步数为:%d\n',kj2);
fprintf('第二种Jacobi迭代cpu时间为:%f\n',tj2);
18.Test44.m
%二维泊松方程，GS迭代实现的两种方式对比
clear,clc;
n=input('请输入二维泊松方程对x,y的分割细密程度（矩阵的维数的平方根）: '); %m为矩阵的维数（二维泊松方程对应的一维的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
T=2*eye(n)+diag(a,1)+diag(a,-1);
I=eye(n);
A=kron(I,T)+kron(T,I);
B=zeros(n,n);
syms x y;
a1(x,y)=x^2+y^2+x; %初值条件
f(x,y)=x^2+x+1+exp(x)+y; %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
Y=linspace(1/(n+1),1-1/(n+1),n); %分割
for i=1:n %构造b
B(i,1)=f(X(i),Y(i))*h^2+a1(X(i),0);
B(1,i)=f(X(i),Y(i))*h^2+a1(0,Y(i));
B(i,n)=f(X(i),Y(i))*h^2+a1(X(i),1);
B(1,n)=f(X(i),Y(i))*h^2+a1(1,Y(i));
end
for i = 2:n-1
for j = 2:n-1
B(i,j)=f(X(i),Y(i))*h^2;
end
end
b=reshape(B,n^2,1);
[tg1,kg1,Ujg1,epsg1]=GS(A,b); %GS:tg是cpu时间,kg是迭代次数,Ujg是A*X=b的解,epsg是误差的范数
t1=1:kg1;
[tg2,kg2,Ujg2,epsg2]=GS1(A,b); %GS1:tg是cpu时间,kg是迭代次数,Ujg是A*X=b的解,epsg是误差的范数
t2=1:kg2;
plot(t1,epsg1,'r*-',t2,epsg2,'b*-')
grid on
xlabel('迭代步数');
ylabel('误差的范数');
title('GS迭代对比图');
legend('GS','GS1')
fprintf('第一种GS迭代步数为:%d\n',kg1);
fprintf('第一种GS迭代cpu时间为:%f\n',tg1);
fprintf('第二种GS迭代步数为:%d\n',kg2);
fprintf('第二种GS迭代cpu时间为:%f\n',tg2);
19.Test55.m
%二维泊松方程，SOR迭代实现的两种方式对比
clear,clc;
n=input('请输入二维泊松方程对x,y的分割细密程度（矩阵的维数的平方根）: '); %m为矩阵的维数（二维泊松方程对应的一维的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
T=2*eye(n)+diag(a,1)+diag(a,-1);
I=eye(n);
A=kron(I,T)+kron(T,I);
w = input('请输入松弛变量w: '); %松弛变量w
B=zeros(n,n);
syms x y;
a1(x,y)=x^2+y^2+x; %初值条件
f(x,y)=x^2+x+1+exp(x)+y; %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
Y=linspace(1/(n+1),1-1/(n+1),n); %分割
for i=1:n %构造b
B(i,1)=f(X(i),Y(i))*h^2+a1(X(i),0);
B(1,i)=f(X(i),Y(i))*h^2+a1(0,Y(i));
B(i,n)=f(X(i),Y(i))*h^2+a1(X(i),1);
B(1,n)=f(X(i),Y(i))*h^2+a1(1,Y(i));
end
for i = 2:n-1
for j = 2:n-1
B(i,j)=f(X(i),Y(i))*h^2;
end
end
b=reshape(B,n^2,1);
[ts1,ks1,Ujs1,epss1]=SOR(A,b,w); %SOR:ts是cpu时间,ks是迭代次数,Ujj是A*X=b的解,epss是误差的范数
t1=1:ks1;
[ts2,ks2,Ujs2,epss2]=SOR1(A,b,w); %SOR1:ts是cpu时间,ks是迭代次数,Ujj是A*X=b的解,epss是误差的范数
t2=1:ks2;
plot(t1,epss1,'r*-',t2,epss2,'b*-')
grid on
xlabel('迭代步数');
ylabel('误差的范数');
title('SOR迭代对比图');
legend('SOR','SOR1')
fprintf('第一种SOR迭代步数为:%d\n',ks1);
fprintf('第一种SOR迭代cpu时间为:%f\n',ts1);
fprintf('第二种SOR迭代步数为:%d\n',ks2);
fprintf('第二种SOR迭代cpu时间为:%f\n',ts2);
20.Test66.m
clear,clc;
n=input('请输入二维泊松方程对x,y的分割细密程度（矩阵的维数的平方根）: '); %m为矩阵的维数（二维泊松方程对应的一维的维数），同时(n+1)也是对区间[0,1]的均匀分割数目
h=1/(n+1); %h为分割后每一个小区间的长度
a=-ones(n-1,1); %构造一维泊松方程的系数矩阵，记为A
T=2*eye(n)+diag(a,1)+diag(a,-1);
I=eye(n);
A=kron(I,T)+kron(T,I);
B=zeros(n,n);
syms x y;
a1(x,y)=x^2+y^2+x; %初值条件
f(x,y)=x^2+x+1+exp(x)+y; %这里是可以改变的，本来想通过键盘输入，但是解析式很难从键盘输入
X=linspace(1/(n+1),1-1/(n+1),n); %分割
Y=linspace(1/(n+1),1-1/(n+1),n); %分割
for i=1:n %构造b
B(i,1)=f(X(i),Y(i))*h^2+a1(X(i),0);
B(1,i)=f(X(i),Y(i))*h^2+a1(0,Y(i));
B(i,n)=f(X(i),Y(i))*h^2+a1(X(i),1);
B(1,n)=f(X(i),Y(i))*h^2+a1(1,Y(i));
end
for i = 2:n-1
for j = 2:n-1
B(i,j)=f(X(i),Y(i))*h^2;
end
end
b=reshape(B,n^2,1);
ts1=zeros(1,19);
ks1=zeros(1,19);
Ujs1=zeros(n^2,19);
k=1;
for w = 0.1:0.1:1.9
[ts1(k),ks1(k),Ujs1(:,k),epss1(1:ks1(k),k)]=SOR(A,b,w); %SOR:ts是cpu时间,ks是迭代次数,Ujj是A*X=b的解,epss是误差的范数
t1=1:ks1(k);
epss=epss1(1:ks1(k),k);
plot(t1,epss)
xlabel('迭代步数');
ylabel('误差的范数');
title('SOR迭代');
fprintf('w=%f时SOR迭代步数为:%d\n',w,ks1(k));
fprintf('w=%f时SOR迭代cpu时间为:%f\n',w,ts1(k));
k=k+1;
figure
end
w=0.1:0.1:1.9;
as=polyfit(w,ks1,5);
fprintf('w和迭代次数的关系: ');
p=poly2str(as,'w');
disp(p);
x1=0.1:0.02:1.9;
y1=as(1)*x1.^5+as(2)*x1.^4+as(3)*x1.^3+as(4)*x1.^2+as(5)*x1+as(6);
plot(w,ks1,'k*',x1,y1,'r-');
legend('原始数据','五次拟合')
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- 数值分析 jacobi迭代法求解线性方程组 MATLAB程序实现
- 基于matlab的jacobi(雅可比)迭代法求解线性方程组
- Jacobi迭代与SOR迭代求解希尔伯特矩阵
- 基于MATLAB的高等数学 求在（x0,y0）处偏导数 数值
- 数值分析——Newton迭代法求解方程附Matlab程序
- MATLAB 求解符号表达式数值的方法：subs函数
- 遗传算法求解带非线性约束的单目标问题，matlab代码，基于K Deb的论文
- 基于Matlab求解diophantine方程
- 编程语言_matlab自定义函数与代数方程求解
- Matlab数值求解超越方程的根
- 基于Matlab的自动控制原理 微分方程求解
- Jacobi迭代、Gauss_Seidel迭代和最佳因子SOR迭代的比较
- 基于MATLAB的有限元法求解EIT技术的正问题
- Jacobi迭代算法matlab版
- Matlab求解代数多项式方程组
- JacobiAndGauss-Seidel迭代(基于java)
- 数值计算-线性方程组求解（2）-追赶法解三对角矩阵-MATLAB实现
- MATLAB线性方程组的迭代求解法
- 基于matlab的Guass-Seidel(高斯--赛德尔) 迭代法求解线性方程组
- 数值计算-线性方程组求解（1）-LU分解-MATLAB实现