//************************************************
//************************************************
// Numerical Algorithms for Visual Computing II
// Sample Solutions for Assignment 2
//
// Copyright Kai Hagenburg 2008
// for code adjustments and bugreporting:
// hagenburg@mia.uni-saarland.de
//************************************************
//************************************************




//************************************************
// Compute the initial vector for exercise 1
// u is the output vector with size 11, i.e. 
// has values u(1)...u(11)
//************************************************

function u = getVector()
  for i=1:6
    u(i)= 2*(i-1)+1
end
  for i=7:11
     u(i) = -7/5*(i-1)+18
end
endfunction


//************************************************
// Additional functions for exercise No. 1
// 
// explicit_time_marching: takes a vector u, a given
// timestep size tau, and a certain number of iterations iter
// sty serves as a variable for different color outputs
// Also the function gives out the result u (Part a)

// Jacobi_iteration: takes a vector u, and a certain
// number of iterations iter.
// sty serves as a variable for different color outputs
// Also the function gives out the result u (Part b)
//************************************************
function u=explicit_time_marching(u,tau,iter,sty)
  // compute length of the vector u
  n = length(u)
  
  // We could hand over delta_x as well, but we 
  // have already predetermined this as h=delta x = 1
  delta_x = 1
  
  // Plot the vector u on the x-axis from 0-10 with style=1
  plot2d(0:10,u,style=1)
  
  // do the iteration over all vector entries
  // also make sure that the fixed boundary conditions are included
  for i=1:iter
    tmp(1) = 1
    tmp(11) = 4
    for j=2:n-1
      tmp(j)= (1-2*tau/delta_x^2)*u(j) + tau/delta_x^2 *u(j-1) + tau/delta_x^2*u(j+1)
    end
    u=tmp
  end
  // Plot the output vector u on the x-axis from 0-10 with style=sty
  plot2d(0:10,u,style=sty)
endfunction

function u=Jacobi_iteration(u,iter,sty)
  // compute length of the vector u
  n=length(u)
  
  // Plot the vector u on the x-axis from 0-10 with style=1  
  plot2d(0:10,u,style=1)
  
  // do the iteration over all vector entries
  // also make sure that the fixed boundary conditions are included  
  for i=1:iter
     tmp(1)=1
     tmp(n)=4
     for j=2:n-1
       tmp(j)=(u(j-1)+u(j+1))/2
     end
     u=tmp
  end
  // Plot the output vector u on the x-axis from 0-10 with style=sty
   plot2d(0:10,u,style=sty)
 endfunction
 
//******************************************************
// This automatically starts the computation for the
// solution of exercise 1
// Note that we used in our computation for the 
// explicit time marching algorithm the use of 
// tau = delta t = 0.25 as seen in the lecture.
//******************************************************
 
function exercise1()
  u=getVector()
  xset("window",1)
  explicit_time_marching(u,0.25,1,2)
  explicit_time_marching(u,0.25,10,3)
  explicit_time_marching(u,0.25,20,4)
  xset("window",2)
  Jacobi_iteration(u,1,2)
  Jacobi_iteration(u,10,3)
  Jacobi_iteration(u,20,4)
endfunction


//*************************************************
// Exercise 2
// buildmatrix builds matrix A, and solution vector u
// according to our calculations for our problem
// exercise2() gives te direct solution for the problem
// by use of the Thomas algorithm
//*************************************************



function [u,A]=buildmatrix()
v(1:9)=-2
A=diag(v)
A(1,2)=1
A(9,8)=1
for i=2:8
  A(i,i-1)=1
  A(i,i+1)=1
end

u(1:9)= 0
u(1)=-1
u(9)=-4
endfunction

function x=exercise2()
[d,A]=buildmatrix()

// simple variant
// v = A\u

// complex variant
// Thomas Algorithm

// Compute length of solution vector
n = length(d)


// Create help vectors a,b,c for the diagonals, i.e.
// a u_i-1 + b u_i + c u_i+1 = d_i
    b(1)=A(1,1)
    b(n)=A(n,n)
    c(1)=A(1,2)
    c(n)=0
    a(n)=A(n,n-1)
    
  for j=2:n-1
    b(j)=A(j,j)
    a(j)=A(j,j-1)
    c(j)=A(j,j+1)
end


// Forward sweep

c(1)=c(1)/b(1)
d(1)=d(1)/b(1)
for i=2:n
  id = (b(i)-c(i-1)*a(i))
  c(i)=c(i)/id
  d(i)=(d(i)-d(i-1)*a(i))/id
end

// Backward substitution

x(n)=d(n)
for i=n-1:-1:1
  x(i)=d(i)-c(i)*x(i+1)
end
  
endfunction
    
    
//****************************************************
// Exercise 3
// needs timestep size tau and number of iterations iter
//****************************************************
    
    

function x=exercise3(tau, iter)

// prepare start vector
v=getVector()
x=v(2:10)
a=v(1)
b=v(11)

// compute system matrix
for i=1:9
  A(i,i)=1+2*tau
  if (i>1) A(i,i-1)=-tau
    end
  if (i<9) A(i,i+1)=-tau
    end
end


// iteratively solve the system
for i=1:iter
  x(1)=x(1)+tau*a
  x(9)=x(9)+tau*b
  x=A\x
end

// plot the result
out(1)=a
out(11)=b
out(2:10)=x
plot2d(1:11,out,style=2)
endfunction


  
