% Harj. 7 ratkaisuja
% 
% teht. 1
%%%%%%%%%%
% m : sisäsolmujen lkm. välin [a,b] jaossa, askel h=1/(m+1)

help spdiags

...

Some elements of B, corresponding to positions "outside" of A, are
    not actually used.  They are not referenced when B is an input and
    are set to zero when B is an output.  If a column of B is longer than
    the diagonal it's representing, elements of super-diagonals of A
    correspond to the lower part of the column of B, while elements of
    sub-diagonals of A correspond to the upper part of the column of B.

% Kokeillaan spdiags:n toimintaa:

yks2345=(1:5)';
A=spdiags([yks2345,yks2345,yks2345],-1:1,5,5);
full(A)
ans =
 
     1     2     0     0     0
     1     2     3     0     0
     0     2     3     4     0
     0     0     3     4     5
     0     0     0     4     5

% Nähdään, että spdiags ottaa aladiagonaalit alusta ja
% ylädiagonaalit lopusta.
% Jos nauha on vakio, tämä käytös ei näy:


%%%%%%%% Ensin blokkityyli, kron:n avulla %%%%%
m=3;
e=ones(m,1);
D=spdiags([e,-2*e,e],-1:1,m,m);
full(D)

ans =
    -2     1     0
     1    -2     1
     0     1    -2

 I=eye(m,m);
 
 
>> help kron

 KRON   Kronecker tensor product.
    KRON(X,Y) is the Kronecker tensor product of X and Y.
    The result is a large matrix formed by taking all possible
    products between the elements of X and those of Y.   For
    example, if X is 2 by 3, then KRON(X,Y) is
 
       [ X(1,1)*Y  X(1,2)*Y  X(1,3)*Y
         X(2,1)*Y  X(2,2)*Y  X(2,3)*Y ]
 
    If either X or Y is sparse, only nonzero elements are multiplied
    in the computation, and the result is sparse.
    
full(kron(D,I))
ans =
    -2     0     0     1     0     0     0     0     0
     0    -2     0     0     1     0     0     0     0
     0     0    -2     0     0     1     0     0     0
     1     0     0    -2     0     0     1     0     0
     0     1     0     0    -2     0     0     1     0
     0     0     1     0     0    -2     0     0     1
     0     0     0     1     0     0    -2     0     0
     0     0     0     0     1     0     0    -2     0
     0     0     0     0     0     1     0     0    -2    
     
full(kron(I,D))     
ans =
    -2     1     0     0     0     0     0     0     0
     1    -2     1     0     0     0     0     0     0
     0     1    -2     0     0     0     0     0     0
     0     0     0    -2     1     0     0     0     0
     0     0     0     1    -2     1     0     0     0
     0     0     0     0     1    -2     0     0     0
     0     0     0     0     0     0    -2     1     0
     0     0     0     0     0     0     1    -2     1
     0     0     0     0     0     0     0     1    -2

%%%%
m=3;
e=ones(m,1);
D=spdiags([e,-2*e,e],-1:1,m,m);
I=eye(m,m);
A=kron(I,D)+kron(D,I);
full(A)
ans =
 
    -4     1     0     1     0     0     0     0     0
     1    -4     1     0     1     0     0     0     0
     0     1    -4     0     0     1     0     0     0
     1     0     0    -4     1     0     1     0     0
     0     1     0     1    -4     1     0     1     0
     0     0     1     0     1    -4     0     0     1
     0     0     0     1     0     0    -4     1     0
     0     0     0     0     1     0     1    -4     1
     0     0     0     0     0     1     0     1    -4
%% Kaunista %%%%


%%% Toinen tapa: suoraan diagonaaleittain

m=3;
M=m^2;
e=ones(M,1);
e1=ones(M-1,1);
mod(1:M-1,m)
e1=mod(1:M-1,m)~=0
e1=e1';
   % Otetaan selvyyden vuoksi täyttöalkioksi jokin "älytön".
A=spdiags([[e1;-10],-4*e,[-10;e1]],-1:1,M,M);
   % Lisätään vielä sivuykköset.
A=A+spdiags([e,e],[-m,m],M,M);
>> full(A)
 
ans =
 
    -4     1     0     1     0     0     0     0     0
     1    -4     1     0     1     0     0     0     0
     0     1    -4     0     0     1     0     0     0
     1     0     0    -4     1     0     1     0     0
     0     1     0     1    -4     1     0     1     0
     0     0     1     0     1    -4     0     0     1
     0     0     0     1     0     0    -4     1     0
     0     0     0     0     1     0     1    -4     1
     0     0     0     0     0     1     0     1    -4

%% Kaunista on tämäkin!  %%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Teht. 2, Moler PDE 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  u'' = exp(-x^2) , x \in [-1,1]

% Numeerisesti

clf
n = 63;         % sisäsolmujen lkm.
h = 2/(n+1);    % Osavälejä n+1
x = (-1:h:1)';
f = exp(-x.^2);
e = ones(n,1);
u = zeros(n+2,1);   % Reunat mukana
A = spdiags([e -2*e e],[-1 0 1],n,n)/h^2;
sisa_i = 2:n+1;
u(sisa_i) = A\f(sisa_i);    % Loksautetaan u:n "sisuksiin" (reunojen sisään)
plot(x,u);
residuaali = max(abs(diff(u,2)/h^2 - exp(-x(sisa_i).^2)))

% Symbolisesti

maple('dsolve({diff(u(x),x,x)=exp(-x^2),u(-1)=0,u(1)=0},u(x));')

ans =

u(x) = 1/2*pi^(1/2)*x*erf(x)+1/2*exp(-x^2)-1/2*erf(1)*pi^(1/2)-1/2*exp(-1)

% Tässä hiukan editointia (pisteet: .* ja .^ )
U=1/2*pi^(1/2)*x.*erf(x)+1/2*exp(-x.^2)-1/2*erf(1)*pi^(1/2)-1/2*exp(-1);


% Tarkan ja numeerisen erotuksen max hilapisteissä:

maksvirhe = max(abs(U - u))

5.1453e-005


%%%%%%%%%%%%%%%%%%%%%
% Teht. 3
%%%%%%%%%%%

nelio=[-3 -3 6]

[u,X,Y]=Poisson0(nelio,@peaks);

contourf(X,Y,u);

subplot(1,2,1)

F=peaks(X,Y);
a=min(min(F)); b=max(max(F));
contourf(X,Y,F,a:1:b)
caxis([a b])
axis square
set(gca,'ytick',get(gca,'xtick'))

subplot(1,2,2)
a = min(min(u));
b = max(max(u));
c = a + (b-a)*(0:.1:10)';
contourf(X,Y,u,c)
caxis([a b])
 axis square
colorbar
set(gca,'ytick',get(gca,'xtick'))

figure
subplot(1,2,1)
meshc(X,Y,F)

subplot(1,2,2)
surfc(X,Y,u)
colorbar

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Teht. 4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Kokeillaan geometrioita:
clear; close all
n=5; z=exp((0:n)*i*2*pi/n);
plot(z)

xv=real(z); yv=imag(z);
fill(xv,yv,'y')
N=9;
N=2*N
h=2/N;
[X,Y]=meshgrid(-1:h:1);
[in,on]=inregion(X,Y,xv,yv);
p = find(in-on);
q = find(on);
hold on
plot(xv,yv,'-',X(p),Y(p),'ko',X(q),Y(q),'ro')
figure
spy(in)
hold on
spy(on)




    

