test

require 'mathpack'

module Mathpack
  module SLE
    # метод прогонки, все коэффиценты прогонки с положительными знаками
    def self.simple_sweep(params = {})
      a = params[:left].dup
      b = params[:right].dup 
      c = params[:diagonal].dup
      f = params[:f].dup
      n = f.length
      fail "left and right should have dimension #{n - 1}" if a.length != n - 1 || b.length != n - 1
      fail "diagonal and f should have dimension #{n}" if c.length != n || f.length != n
      alpha = Array.new(n, b[0] / c[0].to_f)
      beta = Array.new(n, f[0] / c[0].to_f)
      1.upto(n - 2) do |k|
        del = c[k] - alpha[k - 1] * a[k - 1]
        raise 'Division by zero' if del.zero?
        alpha[k] = b[k] / del   
        beta[k] = (f[k] - beta[k - 1] * a[k - 1]) / del
      end
      result = Array.new(n, (f[n - 1] - beta[n - 2] * a[n - 2]) / (c[n - 1] - alpha[n - 2] * a[n - 2]))
      (n - 2).downto(0) do |i|
        result[i] = beta[i] - alpha[i] * result[i + 1]
      end
      result
    end
  end
end

# метод средних прямоугольников
# применим для систем ДУ любой размерности
def explicit_rectangle(params = {})
  nodes = Mathpack::Approximation.generate_nodes(from: params[:from], to: params[:to], step: params[:step])
  y = Array.new(params[:system].length) { Array.new(nodes.length) }
  y.each_index do |i|
      y[i][0] = params[:y0][i]
  end
  for i in 0...nodes.length-1 do
    column = y.transpose[i]
    y_volna = Array.new(y.length)
    for j in 0...y_volna.length do
      y_volna[j] = column[j] + 0.5 * params[:step] * params[:system][j].call(nodes[i], *column)
    end
    for j in 0...y.length do
      y[j][i + 1] = y[j][i] + params[:step] * params[:system][j].call(nodes[i] + 0.5 * params[:step], *y_volna)
    end  
  end
  y
end

# вычисление значения пси
def psi(p, q, r, u1, u2)
  p * u1 + q * u2 - r
end

# метод стрельбы
def shooting_method(p, q, r, left, right, start_condition, cauchie_problem, h)
  n0 = 0.0
  # решаем систему с начальным значением n0
  result1 = explicit_rectangle(from: left, to: right, y0: [n0, start_condition], system: cauchie_problem, step: h)
  # вычислем пси
  psi0 = psi(p, q, r, result1[0].last, result1[1].last)
  # повторяем то же самое с начальным условием n1
  n1 = 5.0
  result2 = explicit_rectangle(from: left, to: right, y0: [n1, start_condition], system: cauchie_problem, step: h)
  psi1 = psi(p, q, r, result2[0].last, result2[1].last)
  # находим истинное значение начального условия и решаем задачу с ним
  n2 = n0 - (n1 - n0) * psi0 / (psi1 - psi0)
  return explicit_rectangle(from: left, to: right, y0: [n2, start_condition], system: cauchie_problem, step: h).first
end

def K(x)
  (x + 2.0)/(x + 3.0)
end

def rs_method(left, right, h)
  a = [] # массивы для коэффициентов прогоночной матрицы
  b = []
  c = []
  f = []
  # вычисляем значения узлов
  nodes = Mathpack::Approximation.generate_nodes(from: left, to: right, step: h)
  n = nodes.length - 1
  # находим значения прогоночных коэффициентов
  b[0] = K(left+0.5*h)/h
  c[0] = -K(left+0.5*h)/h + 0.5 * h * (Math.sin(1.0) - 1.0)*(1.0 + 0.25*h)
  f[0] = -h*(1+0.25*h)
  for i in 1...n do
    a[i-1] = 1.0/(h**2) * (nodes[i] - 0.5*h + 2.0)/(nodes[i] - 0.5*h + 3.0)
    b[i]   = 1.0/(h**2) * (nodes[i] + 0.5*h + 2.0)/(nodes[i] + 0.5*h + 3.0)
    c[i]   = -1.0/(h**2) * (nodes[i] + 0.5*h + 2.0)/(nodes[i] + 0.5*h + 3.0) - 1.0/(h**2) * (nodes[i] - 0.5*h + 2.0)/(nodes[i] - 0.5*h + 3.0) - (1.0 + Math.sin(nodes[i]))
    f[i]   = nodes[i] - 1.0
  end

  c[n]   = K(right-0.5*h)/h + 3.0/8.0 + 0.5*h*(1.0 + Math.sin(1.0))
  f[n]   = 1.5
  a[n-1] = -K(right-0.5*h)/h

  # находим решение системы методом прогонки
  return Mathpack::SLE.simple_sweep(left: a, right: b, diagonal: c, f: f)
end

left = -1.0
right = 1.0
h = 0.1
p = 0.5
q = 1.0
r = 2.0

# система уравнений ДУ для решения методом прогонки
boundary_problem = [ ->(x, u1, u2) { u2 }, -> (x, u1, u2) { (1.0 + Math.sin(x)) * (3.0 + x) / (2.0 + x) * u1 - 1.0 / ((2.0 + x) * (3.0 + x)) * u2 + (x - 1.0) * (x + 3.0) / (x + 2.0) }]

shooting_result = shooting_method(p, q, r, left, right, 0.0, boundary_problem, h)
shooting_result_2 = shooting_method(p, q, r, left, right, 0.0, boundary_problem, 2*h)
Mathpack::IO.print_table_function(filename: 'shooting.csv', x: Mathpack::Approximation.generate_nodes(from: left, to: right, step: h), y: shooting_result)
p Mathpack::IO.count_diff(shooting_result, shooting_result_2) / (2**2 - 1)

rs_result = rs_method(left, right, h)
rs_result_2 = rs_method(left, right, 2*h)
Mathpack::IO.print_table_function(filename: 'rs.csv', x: Mathpack::Approximation.generate_nodes(from: left, to: right, step: h), y: rs_result)
p Mathpack::IO.count_diff(rs_result, rs_result_2) / (2**2 - 1)