Тексты функций

function main(M,N,mu,MAX)

% main(M,N,mu,eps,MAX)

% M - the X-size of the grig

% N - the Y-size of the grig

% MAX - the maximum number of the iterations

% mu - the parameter of the convergence (1<mu<2)

%warning off

%close all

%format short

a = 1/2;

b = 2;

h = a/M;

t = b/N;

q = t/h;

eps=((h+t)/2)^4/1.25;

% initialization

u = ones(M+1, N+1);

x = zeros(M+1,1);

y = zeros(N+1,1);

f = zeros(M+1, N+1);

for i = 1:(M+1),

x(i) = (i-1)*h;

end

for j = 1:(N+1),

y(j) = (j-1)*t;

end

for i=1:M+1

for j=1:N+1

f(i,j)= -(-4*pi^2*y(j)^2+8*pi^2*y(j)+2)*sin(2*pi*x(i));

end;

end;

% the approximation of the boundary conditions

for i = 1:(M+1),

u(i,1) = 0;

u(i,N+1) = 0;

end

for j = 1:(N+1),

u(1,j) = 0;

u(M+1,j) = 0;

end

k = 0;

A = 100;

a1 = 1/(2*(1+q^2));

a2 = q^2;

a3 = t^2*a1;

while (A > eps)

k = k + 1;

A = 0;

for i=2:M,

for j=2:N,

X = a1 * (u(i,j-1) + u(i,j+1) + a2*u(i-1,j) + a2*u(i+1,j)) + a3*f(i,j);

X = mu*X + (1-mu)*u(i,j);

R = abs(X - u(i,j));

if (R > A) A = R; end

u(i,j) = X;

end

end

if k == MAX

sprintf('%d steps were made!!!', MAX)

break;

end

end

sprintf('Quality of approximation is %d', A)

if (k ~= MAX)

sprintf('Riched after %d steps', k)

end

%вывод результатов

y=zeros(M+1,11);

for i=1:M+1

for j=1:11

y(i,j)=u(i,(j-1)*N/10+1);

end;

end;

i=0:1:M;

plot(i/M,y);

legend('t=0','t=1/10 T','t=2/10 T','t=3/10 T','t=4/10 T','t=5/10 T','t=6/10 T','t=7/10 T','t=8/10 T','t=9/10 T','t=T','Location','NorthEastOutside')

title('Параметрический график при различных y');

xlabel('Ось X ');

ylabel('Ось U ');

pause;

X = 0:h:a; %массив делений по оси X

Y = 0:t:b; %массив делений по оси T

z = zeros(21,1); %массив высот

z = min(min(u)):((max(max(u))-min(min(u)))/20):max(max(u));

contour(Y,X,u,z);

colorbar;

title('График изолиний');

grid;

xlabel('Ось Y ');

ylabel('Ось X ');

pause;

mesh(u);

colorbar;

title('3D график полученного решения U(x,y)');

xlabel('Ось Y ');

ylabel('Ось X ');

zlabel('Ось U');

m=1;

n=1;

for i=1:M/10:M+1

for j=1:N/5:N+1

result(m,n)=u(i,j);

n=n+1;

end

n=1;

m=m+1;

end

result

function model(M,N,mu,MAX)

warning off

close all

format short

a = 1;

b = 1;

h = a/M;

t = b/N;

q = t/h;

eps=((h+t)/2)^4/1.25;

% initialization

U = ones(M+1, N+1);

Ur = ones(M+1, N+1);

x = zeros(M+1,1);

y = zeros(N+1,1);

% the approximation of the boundary conditions

for i = 1:(M+1),

x(i) = (i-1)*h;

U(i,1) = cos(x(i));

U(i,N+1) = exp(b)*cos(x(i));

end

for j = 1:(N+1),

y(j) = (j-1)*t;

U(1,j) = exp(y(j));

U(M+1,j) = exp(y(j))*cos(a);

end

for i=1:(M+1)

for j=1:(N+1)

f(i,j)=0;

Ur(i,j)=cos(x(i))*exp(y(j));

end

end

k = 0;

A = 100;

a1 = 1/(2*(1+q^2));

a2 = q^2;

a3 = t^2*a1;

while (A > eps)

k = k + 1;

A = 0;

for i=2:M,

for j=2:N,

X = a1 * (U(i,j-1) + U(i,j+1) + a2*U(i-1,j) + a2*U(i+1,j)) + a3*f(i,j);

X = mu*X + (1-mu)*U(i,j);

R = abs(X - U(i,j));

if (R > A) A = R; end

U(i,j) = X;

end

end

if k == MAX

sprintf('%d steps were made!!!', MAX)

break;

end

end

sprintf('Quality of approximation is %f', A)

if (k ~= MAX)

sprintf('Riched after %d steps', k)

end

%вывод результатов

y1=zeros(M+1,6);

y2=zeros(M+1,6);

for i=1:M+1

for j=1:6

y1(i,j)=U(i,(j-1)*N/5+1);

y2(i,j)=Ur(i,(j-1)*N/5+1);

end;

end;

i=0:1:M;

subplot(2,1,1)

plot(i/M,y1);

legend('t=0','t=1/5 T','t=2/5 T','t=3/5 T','t=4/5 T','t=T','Location','NorthEastOutside')

title('Параметрический график при различных t найденной функции');

xlabel('Ось X ');

ylabel('Ось U ');

subplot(2,1,2)

plot(i/M,y2);

legend('t=0','t=1/5 T','t=2/5 T','t=3/5 T','t=4/5 T','t=T','Location','NorthEastOutside')

title('Параметрический график при различных t модельной функции');

xlabel('Ось X ');

ylabel('Ось U ');

pause;

X = 0:h:a; %массив делений по оси X

Y = 0:t:b; %массив делений по оси T

z = zeros(21,1); %массив высот

z = min(min(U)):((max(max(U))-min(min(U)))/20):max(max(U));

subplot(2,1,1)

contour(Y,X,U,z);

colorbar;

title('График изолиний найденной функции');

grid;

xlabel('Ось T ');

ylabel('Ось X ');

subplot(2,1,2)

contour(Y,X,Ur,z);

colorbar;

title('График изолиний модельной функции');

grid;

xlabel('Ось T ');

ylabel('Ось X ');

pause;

subplot(1,1,1)

mesh(U);

colorbar;

title('3D график найденной функции');

xlabel('Ось T ');

ylabel('Ось X ');

zlabel('Ось U');

pause;

mesh(Ur);

colorbar;

title('3D график модельной функции');

xlabel('Ось T ');

ylabel('Ось X ');

zlabel('Ось U');

maxerror=max(max(abs(Ur-U)));

sprintf('Погрешность вычислений равна %d',maxerror)

function progonka(M, N)

a = 1/2;

b = 2;

h = a/M;

t = b/N;

q = (t/h)^2;

x = zeros(M+1,1);

y = zeros(N+1,1);

f = zeros(M+1, N+1);

func = zeros(M+1, N+1);

for i = 1:(M+1),

x(i) = (i-1)*h;

end

for j = 1:(N+1),

y(j) = (j-1)*t;

end

for i=1:M+1

for j=1:N+1

func(i,j)= -(-4*pi^2*y(j)^2+8*pi^2*y(j)+2)*sin(2*pi*x(i));

end;

end;

C=zeros(M-1,N-1);

C(1,1)=2*(1+q);

C(1,2)=-q;

for j=2:N-2

C(j,j)=2*(1+q);

C(j,j-1)=-q;

C(j,j+1)=-q;

end;

C(M-1,N-1)=2*(1+q);

C(M-2,N-1)=-q;

C(M-1,N-2)=-q;

f=zeros(M-1,1,N-1);

for j=1:N-1

for i=1:M-1

f(i,:,j)=t^2*func(i+1,j+1);

end;

end;

E=zeros(M-1,N-1);

for i=1:M-1

E(i,i)=1;

end;

A=zeros(M-1,N-1,N);

R(:,:,1)=E;

for j=1:N-1

A(:,:,j+1)=(C-A(:,:,j))^-1;

end;

B=zeros(M-1,1,N);

%for i=1:M-1,

B(:,:,1)=0;

%end;

for j=1:N-1

B(:,:,j+1)=A(:,:,j+1)*(B(:,:,j)+f(:,:,j));

end;

u=zeros(M-1,1,N);

for i = 1:M-1,

u(i,:,1) = 0;

u(i,:,N) = 0;

end;

for j=N-1:(-1):1

u(:,:,j)=A(:,:,j+1)*u(:,:,j+1)+B(:,:,j+1);

end;

U=zeros(M+1,N+1);

for i=2:M

for j=2:N

U(i,j)=u(i-1,:,j-1);

end;

end;

for i = 1:(M+1),

U(i,1) = 0;

U(i,N+1) = 0;

end

for j = 1:(N+1),

U(1,j) = 0;

U(M+1,j) = 0;

end

%вывод результатов

y=zeros(M+1,11);

for i=1:M+1

for j=1:11

y(i,j)=U(i,(j-1)*N/10+1);

end;

end;

i=0:1:M;

plot(i/M,y);

legend('t=0','t=1/10 T','t=2/10 T','t=3/10 T','t=4/10 T','t=5/10 T','t=6/10 T','t=7/10 T','t=8/10 T','t=9/10 T','t=T','Location','NorthEastOutside')

title('Параметрический график при различных y');

xlabel('Ось X ');

ylabel('Ось U ');

pause;

X = 0:h:a; %массив делений по оси X

Y = 0:t:b; %массив делений по оси T

z = zeros(21,1); %массив высот

z = min(min(U)):((max(max(U))-min(min(U)))/20):max(max(U));

contour(Y,X,U,z);

colorbar;

title('График изолиний');

grid;

xlabel('Ось Y ');

ylabel('Ось X ');

pause;

mesh(U);

colorbar;

title('3D график полученного решения U(x,y)');

xlabel('Ось Y ');

ylabel('Ось X ');

zlabel('Ось U');

m=1;

n=1;

for i=1:M/10:M+1

for j=1:N/5:N+1

result(m,n)=U(i,j);

n=n+1;

end

n=1;

m=m+1;

end

result

function y = epsilon(M,N,MAX,mu)

a = 1/2;

b = 2;

h = a/M;

t = b/N;

q = t/h;

x = zeros(M+1,1);

y = zeros(N+1,1);

u = ones(M+1, N+1);

f = zeros(size(u));

for i = 1:(M+1),

x(i) = (i-1)*h;

end

for j = 1:(N+1),

y(j) = (j-1)*t;

end

for i=1:M+1

for j=1:N+1

f(i,j)=-(-4*pi^2*y(j)^2+8*pi^2*y(j)+2)*sin(2*pi*x(i));

end;

end;

for i = 1:(M+1),

u(i,1) = 0;

u(i,N+1) = 0;

end

for j = 1:(N+1),

u(1,j) = 0;

u(M+1,j) = 0;

end

k = 0;

A = 100;

a1 = 1/(2*(1+q^2));

a2 = q^2;

a3 = t^2*a1;

for k = 1:MAX,

A = 0;

for i=2:M,

for j=2:N,

X = a1*(u(i,j-1)+u(i,j+1)+a2*u(i-1,j)+a2*u(i+1,j))+a3*f(i,j);

X = mu*X + (1-mu)*u(i,j);

R = abs(X - u(i,j));

if (R > A) A = R; end

u(i,j) = X;

end

end

end

y=A;

y=A;