% Lämpöyhtälön numeriikkaa
% 6.3.05 HA   matlab/lampoyhtnumeriikkaa.m

%Differenssimenetelmä

% 1) Eksplisiittinen (huono)
% Voidaan toki kirjoittaa helposti Matlab-funktio (kts. lampodiffE).
% Toisaalta, kun tulimme kirjoittaneeksi Maple-funktion, siirretään tulosmatriisi 
%   ==> saadaan käyttöön Matlabin tehokkaammat visualisointivälineet.
clear ; close all
V=load('C:/MATLAB6p1/work/V.txt');
surfc(V); colorbar; shg
subplot(1,2,1)
plot(V); 
subplot(1,2,2)
plot(V')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cooper s. 127 tehäviä
%   Sovelletaan LYdiffB-funktiotamme
% Teht. 11
% 
f=inline('(x>0 & x<5).*x + (x>=5 & x<=10 ).*(10-x)','x')
g=inline('zeros(size(t))','t');
% g=inline('0','t')  ei ole "array smart"

nx=60; nt=100; dt=0.05;
vali=[0 10];
[x,t,U]=LYdiffB(f,g,g,vali,nx,nt,dt);
figure(1), clf, surfc(x,t,U)
xlabel('x');ylabel('t');zlabel('lämpötila u')


snap1=find(t==1)
plot(x,U(snap1,:))

% Teht. 12
%
clear; close all
f=inline('-sin(pi*x/5).*exp(x/2)','x');
fplot(f,[0 10])
g=inline('zeros(size(t))','t');
nx=60; nt=100; dt=0.05;
nx=120; nt=200; dt=0.05;
vali=[0 10];
[x,t,U]=LYdiffB(f,g,g,vali,nx,nt,dt);
[x,t,V]=CrNich(f,g,g,vali,nx,nt,dt);
ero=max(max(U-V))
figure; meshc(x,t,V)
xlabel('x');ylabel('t');zlabel('lämpötila u')
figure
plot(x,V(nt,:));shg

maxarvot=max(V');
[M,N]=size(V);
maxmat=repmat(maxarvot',1,N);
[X,T]=meshgrid(x,t);
xmax=X(maxmat==V);
tmax=T(maxmat==V);
plot(xmax,tmax)
plot(xmax,tmax,'.')

xlabel('x');ylabel('t')
umax=V(maxmat==V);
plot3(xmax,tmax,umax)
hold on
meshc(x,t,V)
plot3(xmax,tmax,umax,'r')
Uanal=lamposin(1,x,t,'HALb');
plot(x,U);shg

%%%%%%%%%%%%%%%%%%%%%
% Maplella generoitu analyyttinen ratkaisu:
Uanal=lamposin(1,x,t,'HALb'); 

plot(x,U);shg
figure
plot(x,Uanal);shg
plot(x,Uanal-U);shg
W=Uanal-U;

plot(x,W(10:end,:);shg

>> xi=x(1:10:end);
>> ti=t(1:10:end);
>> [Xi,Ti]=meshgrid(xi,ti);
>> [X,T]=meshgrid(x,t);
>> Vi=interp2(X,T,V,Xi,Ti);
>> mesh(Xi,Ti,Vi)
>> shg
>> figure
>> mesh(Xi,Ti,Vi)
>> surfc(Xi,Ti,Vi)
>> surfc(xi,ti,Vi)

     function [Xi,Yi,Zi]=harvenna(x,y,Z,xaskel,yaskel)
     xi=x(1:xaskel:end);
     yi=y(1:yaskel:end);
     [X,Y]=meshgrid(x,y);
     [Xi,Yi]=meshgrid(xi,yi);
     Zi=interp2(X,Y,Z,Xi,Yi);

>> [Xi,Yi,Zi]=harvenna(x,t,V,10,10);
>> surfc(Xi,Yi,Zi)
>> shg



% Teht. 13 s. 128, harj 5 teht. 6

clear; close all; format compact
f=inline('x/2-2+sin(4*pi*x/10)','x');
ga=inline('-2*ones(size(t))','t');
gb=inline('3*ones(size(t))','t');
fplot(f,[0 10])
nx=60; nt=100; dt=0.05; 
%nx=60; nt=100; dt=0.1; T=nt*dt
vali=[0 10];
[x,t,U]=LYdiffB(f,ga,gb,vali,nx,nt,dt);
[x,t,V]=CrNich(f,ga,gb,vali,nx,nt,dt);
subplot(211)
meshc(x,t,U)
xlabel('x');ylabel('t');zlabel('lämpötila u')
subplot(212)
meshc(x,t,V)
xlabel('x');ylabel('t');zlabel('lämpötila u'); shg
figure
plot(x,U(100,:));shg
figure
plot(x,U(90:end,:));shg
figure
plot(x,V(100,:));shg
figure
plot(x,V(90:end,:));shg

%%%%% Mathews:
clear; close all; format compact
f=inline('x/2-2+sin(4*pi*x/10)','x');
ga=inline('-2*ones(size(t))','t');
gb=inline('3*ones(size(t))','t');
fplot(f,[0 10])
nx=60; nt=100; dt=0.05; T=nt*dt
U = Mathcrnich(f,ga,gb,10,nt*dt,1,nx,nt);
%nx=60; nt=100; dt=0.1; T=nt*dt
vali=[0 10];
[x,t,U]=Mathcrnich(f,ga,gb,vali,nx,nt,dt);

% Teht. 14

clear; close all
f=inline('zeros(size(x))','x');
ga=inline('t./(1+t)+t.*exp(-t/5)','t');
gb=inline('zeros(size(t))','t');
nx=60; nt=100; dt=0.05; Tf=nt*dt ; Tsiivut=0:Tf;
%nx=60; nt=100; dt=0.1; Tf=nt*dt
%nx=60; nt=200; dt=0.1; Tf=nt*dt
vali=[0 10];
%[x,t,U]=LYdiffB(f,ga,gb,vali,nx,nt,dt);
[x,t,V]=CrNich(f,ga,gb,vali,nx,nt,dt);
%meshc(x,t,U); title('LYdiffB implisiittinen'); figure 
meshc(x,t,V); title('CrNich')
xlabel('x');ylabel('t');zlabel('lämpötila u')
figure
plot(x,U(nt,:));shg
figure
plot(x,U(nt-10:end,:));shg

figure
hold on
for nt=100:100:500
    [x,t,V]=CrNich(f,ga,gb,vali,nx,nt,dt);
    plot(x,V(nt,:));shg
end
hold off;plot(x,V(nt-20:nt,:),'r');shg

fplot(ga,[0,20])

>> fplot(ga,[30,40]);shg
>> fplot(ga,[40,100]);shg
>> fplot(ga,[100,200]);shg
>> fplot(ga,[200,300]);shg
>> T=linspace(10,40);
>> T=linspace(200,400);
>> y=ga(T);
>> max(y)
ans =
    0.9968
>> y(end)
ans =
    0.9968

    nt=5000; Tf=nt*dt
 [x,t,V]=CrNich(f,ga,gb,vali,nx,nt,dt);
    plot(x,V(nt,:));shg
    
 % Tästäpä tuli kauniisti suora (0,1) -> (10,0)
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Exe 4.3. s. 138 ->
% 20.3.05
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Teht. 7
% Neumannin RE't

% Lämpöyhtälö: c=vakio, m=0, f(x,t,u,dudx)=dudx, s=0
% (AE) : u(x,t0)= u0(x)
function pdelampoNeumann
close all
global x t u     % Ennen kutsua myös istunnossa, jos haluat jälkeenpäin käyttää.
m = 0;
L=10; x = linspace(0,L,20);
%t = linspace(0,1,5);
t = linspace(0,2,10);
sol = pdepe(m,@lampo,@AE,@RE,x,t);
% Extract the first solution component as u.

u = sol(:,:,1);

% A surface plot is often a good way to study a solution.
surf(x,t,u)    
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')
% A solution profile can also be illuminating.
figure
plot(x,u(end,:))
title('Solution at t = 2')
xlabel('Distance x')
ylabel('u(x,2)')
% --------------------------------------------------------------
% Lämpöyhtälö: c=1 (huomaa: c on vasemmalla puolella, eli pi^2 u_t = u_xx
function [c,f,s] = lampo(x,t,u,DuDx)
c = 1;
f = DuDx;
s = 0;
% --------------------------------------------------------------
function u0 = AE(x)
u0 = exp(-(x-5).^2);
% --------------------------------------------------------------
function [pl,ql,pr,qr] = RE(xl,ul,xr,ur,t)
% Neumannin RE:t  
pl = 0;                    % Vasemmassa reunassa derivaatta = 0
ql = 1;     
%pr=0;
pr = 0;                   % Oikeassa reunassa derivaatta u_x = t  
qr = 1;                   
