Menu
Article Homework

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

NodeKey

NodeKey

• 10 min read

内容包括

  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-8
$inf=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-$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*$inf)
      lans=2*k-rans
      puts '测试方程为2x^2+4x+2=18'
      printf("方程解为：%.2f\n",lans)
      printf(" %.2f",rans) if((lans-rans).abs>$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>$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>$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 > $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 > $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 < $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
Share