Title / Description
Code 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)
Author
Highlight as C C++ CSS Clojure Delphi ERb Groovy (beta) HAML HTML JSON Java JavaScript PHP Plain text Python Ruby SQL XML YAML diff code
1.1.2
1.1.2