Ruby

数值计算编程作业(Ruby版)

Posted on

内容包括

  1. 二分法求解二元一次方程的根
  2. 牛顿法求解一元多次方程的根
  3. 顺序消去法解线性方程组
  4. 列选主元消去法解线性方程组
  5. 全选主元消去法解线性方程组
  6. Doolittle分解线性方程组
  7. Crout分解线性方程组
  8. 平方根法求线性方程组
  9. 拉格朗日插值法
  10. 牛顿插值法
  11. 最小二乘法——线性拟合
  12. 变步长梯形求积算法计算积分
  13. 龙贝格算法计算积分
  14. 改进欧拉算法求常微分方程的数值解
  15. 四阶龙格-库塔法求常微分方程的数值解

累!!!!!!!不详细介绍了

#!/usr/bin/env ruby
# encoding: utf-8

eps=1e-8inf=0x3f3f3f3f

class TwoSplit
  def initialize
    @a,@b,@c,@d=2,4,2,18
  end

  def getVal(x)
    @a*x*x+@b*x+@c-@d
  end
  private :getVal

  def three_split
    le,ri,mid=-inf,inf,1
    while le+eps<ri       mid=(le+ri)/2       if(getVal(mid)+eps<0)
        le=mid+eps       else         ri=mid-eps
      end
    end
    mid
  end
  private :three_split

  def two_split(le,ri)
    mid=1
    while le+eps<ri       mid=(le+ri)/2       if(getVal(mid)<0)         le=mid+eps
      else
        ri=mid-
*** QuickLaTeX cannot compile formula:
eps
      end
    end
    mid
  end
  private :two_split

  def calculate_ans
    @a,@b,@c,@d=-@a,-@b,-@c,-@d if (@a<0)
    k=three_split
    if(getVal(k)>0)
      puts 'No Answer'
    else
      rans=two_split(k,1.0*

*** Error message:
Missing $ inserted.

inf)
lans=2*k-rans
puts '测试方程为2x^2+4x+2=18'
printf("方程解为:%.2f\n",lans)
printf(" %.2f",rans) if((lans-rans).abs>

*** QuickLaTeX cannot compile formula:
eps*10)
      puts ''
    end
  end
end

class NewtonIterate
  def initialize
    @a,@b,@c,@d=-2,-5,1,3
  end

  def get_val(x)
    @a*x*x+@b*x+@c-@d
  end

  def get_der(x)
    2*@a*x+@b
  end

  def newton(val)
    pre=val
    val=pre-get_val(pre)/get_der(pre)
    while((pre-val).abs>

*** Error message:
Missing $ inserted.

eps)
pre=val
val=pre-get_val(pre)/get_der(pre)
end
val
end

def calculate_ans
if(@b*@b-4*@a*@c<0)
puts 'No Answer'
else
lans=newton(-1.0*inf)       rans=newton(1.0*inf)
puts '求解方程为-2x^2-5x+1=3'
printf("方程解为:%.2f",lans);
printf(" %.2f",rans) if((lans-rans).abs>

*** QuickLaTeX cannot compile formula:
eps)
      puts ''
    end
  end
end

class ResolveMat
  def initialize
    @mat=6.times.map{ [0] * 6 }
    @l=6.times.map{ [0] * 6 }
    @u=6.times.map{ [0] * 6 }
    @x,@y=[0]*6,[0]*6
    @b=[1,2,3,4]
    @mat[0]=[2,10,4,-3]
    @mat[1]=[-3,-4,-12,13]
    @mat[2]=[1,2,6,-4]
    @mat[3]=[4,14,9,-13]
    @n=4
  end

  def print_mat(mat)
    yield
    @n.times do |i|
      print '['
      @n.times do |j|
        if((mat[i][j]-0).abs >

*** Error message:
Missing $ inserted.

eps)
printf("%8.3f ", mat[i][j]);
else
printf(" ");
end
end
puts ']'
end
puts ''
end

def doolittle
@n.times { |i| @l[i][i]=1 }

@n.times do |k|
k.upto(@n-1) do |j|
@u[k][j]=@mat[k][j]
k.times { |i| @u[k][j]-=(@l[k][i]*@u[i][j]) }
end

(k+1).upto(@n-1) do |i|
@l[i][k]=@mat[i][k]
k.times { |j| @l[i][k]-=(@l[i][j]*@u[j][k]) }
@l[i][k]/=1.0*@u[k][k]
end
end

@n.times do |i|
@y[i]=@b[i]
i.times { |j| @y[i]-=(@l[i][j]*@y[j]) }
end

(@n-1).downto(0) do |i|
@x[i]=@y[i]
(i+1).upto(@n-1) { |j| @x[i]-=(@u[i][j]*@x[j]) }
@x[i]/=1.0*@u[i][i]
end
puts 'Doolittle分解矩阵'
print_mat(@mat) {puts '原矩阵为:'}
print_mat(@l) {puts '矩阵L为:'}
print_mat(@u) {puts '矩阵U为:'}
end

def crout
@n.times { |i| @u[i][i]=1 }

@n.times do |k|
k.upto(@n-1) do |i|
@l[i][k]=@mat[i][k]
k.times { |j| @l[i][k]-=(@l[i][j]*@u[j][k]) }
end
(k+1).upto(@n-1) do |i|
@u[k][i]=@mat[k][i]
k.times { |j| @u[k][i]-=(@l[k][j]*@u[j][i]) }
@u[k][i]/=1.0*@l[k][k]
end
end

@n.times do |i|
@y[i]=@b[i]
i.times { |j| @y[i]-=(@l[i][j]*@y[j]) }
@y[i]/=1.0*@l[i][i]
end
(@n-1).downto(0) do |i|
@x[i]=@y[i]
(i+1).upto(@n-1) { |j| @x[i]-=(@u[i][j]*@x[j]) }
# @x[i]/=@u[i][i]
end
puts 'Crout分解矩阵'
print_mat(@mat) {puts '原矩阵为:'}
print_mat(@l) {puts '矩阵L为:'}
print_mat(@u) {puts '矩阵U为:'}
end

def square #待验证
t=6.times.map { [0]*6 }
lt=6.times.map { [0]*6 }
d=6.times.map { [0]*6 }

@n.times { |i| lt[i][i]=1 }
t[0][0]=@mat[0][0]

1.upto(@n-1) do |i|
t[i][0] = @mat[i][0]
lt[0][i] = @mat[0][i] / t[0][0]
end

1.upto(@n-1) do |i|
1.upto(i) do |j|
t[i][0]=@mat[i][0]
j.times { |k| t[i][j]-=t[i][k]*lt[k][j] }
end

(i+1).upto(@n-1) do |j|
lt[i][j] = @mat[i][j]
i.times {|k| lt[i][j] = lt[i][j] - t[i][k] * lt[k][j] }
lt[i][j] = lt[i][j] / t[i][i]
end
end

@n.times do |i|
@n.times do |j|
@l[i][j]=lt[i][j]
d[i][j]=0 if i!=j
d[i][j]=t[i][j] if i==j
end
end

@y[0]=@b[0]
1.upto(@n-1) do |i|
i.times.reduce(0.0) { |sum,k| @y[i]= @b[i]-(sum+=lt[k][i]*@y[k]) }
# m=0.0
end

@x[@n-1]=@y[@n-1]/t[@n-1][@n-1]
(@n-1).downto(0) do |i|
(i+1).upto(@n-1).reduce(0.0) { |sum,k| @x[i]=@y[i]/d[i][i] - (sum+=lt[i][k]*@x[k]) }
end

print_mat(@mat) { puts '所要求解的线性方程为:'}
@n.times { |i| printf("x[%d]=%.3f \n", i+1, @x[i]); }
end
end

class Order

def initialize
@mat = Array.new(4,[])
@mat[1]=[0,1,2,3,1]
@mat[2]=[0,5,4,10,0]
@mat[3]=[0,3,-0.1,1,2]
@n=3
end

def print_mat
yield
@n.times do | i |
(@n+1).times do | j |
printf("%.2f ",@mat[i+1][j+1])
end
puts ''
end
end
private :print_mat

def print_ans
print "解为:"
@n.times do | k |
printf("%.2f ",@mat[k+1][@n+1])
end
printf "\n\n"
end
private :print_ans

def calculate_answer
@mat[@n][@n+1] /= @mat[@n][@n]
(@n-1).downto(1) do |i|
@mat[i][@n+1]-=(i+1..@n).reduce(0) {|sum,j| sum+=@mat[i][j]*@mat[j][@n+1]}
end
end
private :calculate_answer

def calculate_mat(i)
(i+1).upto(@n+1) do | j |
@mat[i][j]/=(1.0*@mat[i][i])
end
@mat[i][i]=1
(i+1).upto(@n) do | k |
(i+1).upto(@n+1) do | j |
@mat[k][j]-=@mat[k][i]*@mat[i][j]
end
end
end
private :calculate_mat

def original_order
print_mat {puts "顺序消去法:\n原矩阵为:"}
1.upto(@n-1) do | i |
calculate_mat(i)
end
print_mat {puts '所得上三角矩阵为:'}
calculate_answer
print_ans
end

def col_order
print_mat {puts "列主元素消去法:\n原矩阵为:"}
1.upto(@n-1) do | i |
now,maxx=i,@mat[i][i]
(i+1).upto(@n) do | j |
if(@mat[j][i]>maxx)
maxx,now=@mat[j][i],j
end
end

if i!=now
(i..@n+1).each do | k |
@mat[i][k],@mat[now][k] = @mat[now][k],@mat[i][k]
end
end
calculate_mat(i)
end
print_mat {puts '所得上三角矩阵为:'}
calculate_answer
print_ans
end

def all_order
print_mat {puts "全主元素消去法:\n原矩阵为:"}
1.upto(@n-1) do | i |
nx,ny,maxx=i,i,@mat[i][i]

i.upto(@n) do | j |
i.upto(@n) do | t |
if(@mat[j][t]>maxx)
maxx,nx,ny=@mat[j][t],j,t
end
end
end

if i!=nx #row exchange
i.upto(@n+1) do | k |
@mat[i][k],@mat[nx][k] = @mat[nx][k],@mat[i][k]
end
end

if i!=ny # col exchange
i.upto(@n) do | k |
@mat[k][i],@mat[k][ny] = @mat[k][ny],@mat[k][i]
end
end

calculate_mat(i)
end
print_mat {puts '所得上三角矩阵为:'}
calculate_answer
print_ans
end

end

class VarStep
def initialize
@val = (1 + getVal(1))/2
@pre,@n = 0x3f3f3f3f,1
end

def getVal(x)
Math.sin(x) / x
end

def func(n,pre)
sum,x = 0.0,0
n.times do
x = x + 1.0/(2*n)
sum = sum + getVal(x)
x = x + 1.0/(2*n)
end
@pre/2.0 + sum/(2*n)
end

def calculate_ans
while (@val-@pre).abs >

*** QuickLaTeX cannot compile formula:
eps
      @pre = @val
      @val = func(@n,@val)
      @n*=2
    end
    puts '积分函数为sin(x)/x,积分区间为[0,1]'
    printf("积分结果为%.8f\n",@val)
  end

end

class FourthRomberg
  def initialize
    @a,@b,@n=1,5,5
    @x,@y,@k=[],[],[]
    @h=1.0*(@b-@a)/@n
    @x[0],@y[0]=@a,17
  end

  def der_func(x,y)
    x-y+1
  end
  private :der_func

  def calculate_ans
    1.upto(@n) do |i|
      @x[i-1]=@x[0]+(i-1)*@h
      h1=@h/2.0
      @k[1]=der_func(@x[i-1],@y[i-1])
      @k[2]=der_func(@x[i-1]+h1,@y[i-1]+h1*@k[1])
      @k[3]=der_func(@x[i-1]+h1,@y[i-1]+h1*@k[2])
      @k[4]=der_func(@x[i-1]+@h,@y[i-1]+@h*@k[3])
      @y[i]=@y[i-1]+@h*(@k[1]+2*@k[2]+2*@k[3]+@k[4])/6.0
    end

    puts "a=1,b=5\nn=5\ny0=17\nf(x)=x-y+1\n"
    @n.times do |i|
      printf("y[%d]: %.2f\n",i+1,@y[i+1])
    end
  end
end

class ImproveEuler
  def initialize
    @up,@down,@n = 0.5,0,5
    @h = (@up-@down) / @n
    @xn = @down
    @yn = 1
  end

  def der_func(x,y)
    x - y + 1
  end

  def calculate_ans
    print "f(x,y)=x-y+1\n积分区间为[0,0.5]\ny0=1\n区间5等分\n"
    @n.times do | i |
      xl = @xn + @h
      ylb = @yn + @h * der_func(@xn,@yn)
      yl = @yn + @h / 2 * (der_func(@xn,@yn)+der_func(xl,ylb))
      printf("x%d=%f,y%d=%f\n", i, xl, i, yl)
      @xn,@yn = xl,yl
    end
  end
end

class Lagrance
  def initialize
    @n,@m,@ans=2,195,0
    @x=[121,144,100]
    @y=[11,12,10]
  end

  def calculate_ans
    0.upto(@n) do |i|
      key=1.0*@y[i]
      0.upto(@n) do |j|
        next if(i==j)
        key*=(1.0*(@m-@x[j])/(@x[i]-@x[j]))
      end
      @ans+=key
    end
    puts '给定3个坐标:'
    puts 'x[0]:121,y[0]:11'
    puts 'x[1]:144,y[1]:12'
    puts 'x[2]:100,y[2]:10'
    printf("询问195的y值为%.2f\n",@ans)
  end
end

class NewtonInterpolate
  def initialize
    @n=4
    @form=[0.0]*6
    @x=[10,9,3,11,2]
    @y=[100,80,9,123,4]
  end

  def calculate_ans
    @n.times { |i| @form[i]=@y[i] }
    (@n+1).times do |i|
      @n.downto(i+1) { |j| @form[j] = (@form[j]-@form[j - 1]) / (@x[j] - @x[j - 1 - i]) }
    end
    tmp,hx=1,6.42
    newton=@form[0]
    @n.times do |i|
      tmp*=(hx-@x[i])
      newton+=tmp*@form[i+1]
    end
    @n.times { |i| printf("第%d个点%.2f %.2f\n",i+1,@x[i],@y[i]) }
    printf("F(%.2f) = %.2f\n", hx, newton)
  end
end

class LeastSquare
  def initialize
    @x=[10,9,5,6,2]
    @y=[99,81,24,37,4]
    @a,@b=0.0,0.0
    @n=5
  end

  def get_val(x)
    @a*x+@b
  end

  def calculate_ans
    t=[0.0]*4
    @n.times do |i|
      t[0]+=@x[i]*@x[i]
      t[1]+=@x[i]
      t[2]+=@x[i]*@y[i]
      t[3]+=@y[i]
    end
    @a=(t[2]*@n-t[1]*t[3])/(t[0]*@n-t[1]*t[1])
    @b=(t[0]*t[3]-t[1]*t[2])/(t[0]*@n-t[1]*t[1])
    printf("所得函数为%dx%d\n",@a,@b)
    puts '测试数:12'
    puts get_val(12)
  end
end

class Romberg
  @@pi=4.0*Math.atan(1.0)
  def initialize
    @a,@b=0.000001,1.0
    @n,@h=20,Float(@b-@a)
    @t=20.times.map{ [0] * 2 }
    @t[0][1]=trapezium
    @n*=2
  end

  def get_val(x)
    Math.sin(x)/x
  end
  private :get_val

  def trapezium
    @h=(@b-@a)/@n
    sum = 1.upto(@n-1).reduce(0.0) {|tmp,i| tmp+=get_val(@a+i*@h) }
    sum+=(get_val(@a)+get_val(@b))/2.0
    @h*sum
  end

  def calculate_ans

    puts 'sin(x)/x在[0,1]区间内的积分测试:'
    (1..9).each do |m|
      m.times { |i| @t[i][0]=@t[i][1] }
      @t[0][1]=trapezium
      @n*=2
      m.times do |i|
        @t[i+1][1]=@t[i][1]+(@t[i][1]-@t[i][0]) / ((4**m)-1)
      end
      if (@t[m-1][1]-@t[m][1]).abs <

*** Error message:
Missing $ inserted.

eps
printf("计算的数为:%.6f\n", @t[m][1])
return @t[m][1]
end
end
puts 'No Answer'
end
end

def nodekey
while true
puts '<===================================>'
puts '1.二分法求解二元一次方程的根'
puts '2.牛顿法求解一元多次方程的根'
puts '3.顺序去法解线性方程组'
puts '4.列选主元消去法解线性方程组'
puts '5.全选主元消去法解线性方程组'
puts '6.Doolittle分解线性方程组'
puts '7.Crout分解线性方程组'
puts '8.平方根法求线性方程组'
puts '9.拉格朗日插值法'
puts '10.牛顿插值法'
puts '11.最小二乘法——线性拟合'
puts '12.变步长梯形求积算法计算积分'
puts '13.龙贝格算法计算积分'
puts '14.改进欧拉算法求常微分方程的数值解'
puts '15.四阶龙格-库塔法求常微分方程的数值解'
puts ' 0.退出'
puts '<====================================>'
print "请选择你要选择的操作:"
op=gets.to_i
case op
when 0
break
when 1
TwoSplit.new.calculate_ans
when 2
NewtonIterate.new.calculate_ans
when 3
Order.new.original_order
when 4
Order.new.col_order
when 5
Order.new.all_order
when 6
ResolveMat.new.doolittle
when 7
ResolveMat.new.crout
when 8
ResolveMat.new.square
when 9
Lagrance.new.calculate_ans
when 10
NewtonInterpolate.new.calculate_ans
when 11
LeastSquare.new.calculate_ans
when 12
VarStep.new.calculate_ans
when 13
Romberg.new.calculate_ans
when 14
ImproveEuler.new.calculate_ans
when 15
FourthRomberg.new.calculate_ans
else
puts '请输入0-15的选项,选择0为退出'
end
end
end

nodekey