#control problem parabolic i.e.
#we want to find s(t) (approximation of s(t)) such that it minimizes
#\int_0^T\int_{T_1}^T |u(t,x)-h(t,x)|^2dx dt
#where h given, 0<=T_1<T
#u the solution of 
#u_t=u_{xx} 
#u(0,x)=u0(x) - given initial cond
#u(t,0)=s(t) #control
#u'(t,1)=0 zero N bnd
#We discretize the parabolic pdes by FDM on equidistant mesh
#and assume that s is either a polynomial in t or 
#a linear spline (more reasonable)
#we can also change the right bnd to Robin or Neuman bnd i.e .e.g. u'(0,1)=s(t)
#the FDM scheme 
# d/dt u=(u_{k-1}-2u_k + u_{k+1})/h^2 k=1,...,N-1
#u_0=s(t)
#(u_N-u_{N-1})/h=0 --> u_n=u_{N-1}
#taking into account this we get 
#ODE: u'=(-1/h^2)Au +(1/h^2)s(t)e1 (e1 first unit vector)

1==1;


function L=lap1dmix(n)
#creates sparse discrete minus 1d laplacian on
#uniform grid scaled by h^2 i.e. 3 diag matrix with 
#2 on the main diag except L(n,n)=1, 
#-1 on 1st super- and sub- diagonals
#n size of L nxn
#OUTPUT 
#L 3diag matrix
#
#Example: L=lap1dmix(4); full(L) 
 S=-ones(1,3*n-2);
 S(1:n-1)*=-2;S(n)=1;
 I=[1:n,1:n-1,2:n];#row indices: diagonal; superdiag and subdiag
 J=[1:n,2:n,1:n-1];#column indices: diagonal; superdiag and subdiag
 L=sparse(I,J,S,n,n);
end

function y=vf(x,t)
#vector field for ODE: y'=(-1/h^2)*Ax+(1/h^2)[s(t);0,...;0]
#h mesh parameter 
#A given matrix s function handle t oa scalar function
#h,A,ss (function handle to s()) -parametres passed as global vars
#OUTPUT
#y=(-1/h^2)*A*x+(1/h^2)[s(t);0,...;0]
 global h A ss
 y=(-1/(h*h))*(A*x);
 y(1)=y(1)+ss(t)/(h*h);
endfunction


function [u,t,x]=mixparabsolver(u0,s,N=10,M=10,a=0,b=2,T=2)
#solver for mixed parabolic  problem u_t-u''=0 
#u(t,a)=s u'(t,b)=0 by FDM on equidistant mesh
#u_0=s(t) Dirichlet bnd
# u'=(u_{k-1}-2u_k +u_{k+1})/(h*h) k=1..N-1
#u_N-u_{N-1}=h*0 Neumann bnd
#INPUT
#u0 - function handle to initial value
#s - function handle to Dirichlet bnd at a 
#a,b, - interval ends
#N - no of subintervals - N+1 no of meshpoints 
#M - no of time meshpoints (M+1)
#OUTPUT
#u- matrix (M+1) x (N+1)
#   in (k-th) row i.e. u(k,:) the aproximation 
#    of the FDM solution for t(k)=(k-1)*(T/M)
#t - vector of times ((k-1)*T/M) k=0,..,M
#x - vector of space mesh (x(k)=(k-1)*h) k=1,..,N+1; h=(b-a)/N 
x=linspace(a,b,N+1); x=x';
U0=u0(x(2:N)); #initial value
global h A ss
ss=s;
A=sparse(N-1,N-1);
A=lap1dmix(N-1,N-1);
h=(b-a)/N;
t=linspace(0,T,M+1);
u=lsode(@vf,U0,t);
u=[s(t)',u,u(:,N-1)];
endfunction


#test
L=2;
T=1;
n=50;
M=20;N=2*n;
s=@(t) 1-t;
u0=@(x) (x<1);
[u,t,x]=mixparabsolver(u0,s,N,M,a=0,b=L,T);plot(u')
#indx indt indices of proper space points/times
#for which compute difference of u and h
indx=(n+1):2*n;
indt=1:M+1;
xh=x(indx);
th=t(indt);
uh=u(indt,indx);
h=diag(cos(0.5*pi*th))*ones(size(uh));
er=uh-h;

#it remains to write objective function i.e.
#which we want t ominimize i.e. for given 
#u - FDM approx of solution on equidistant meshes 
#t -  time mesh, x space mesh
#er - error on this mesh
#compute approx of double integral using trapz (twice? in a loop?)
 


