連立方程式をJuliaで

本日は

ピポット選択つきのガウス消去法を用いた連立方程式を解いてみるのJulia版です．

C言語による数値計算入門

で勉強していたPythonのコードをJuliaに書き換えたものになります．

実装例(Julia)

#solve linear equation

#Example1
a=float([0 4 5 2
         1 0 2 -6;
         4 1 0 -2;
         1 7 1 0])::Array{Float64,2}
b=float([9,-3, 1, -3])::Array{Float64,1}
x=float([0, 0, 0, 0])::Array{Float64,1}

function main(a,b,x)
    ori_a=a
    a=hcat(a,b)
    row,col=size(a)
    for j in 1:row
        #search pivot
        maxidx=j-1+indmax(abs.(a[j:end,j]))
        #swap
        a[j,:],a[maxidx,:]=a[maxidx,:],a[j,:]
        #push foward
        for i in j+1:row
            p=-a[i,j]/a[j,j]
            a[i,:]+=p*a[j,:]
        end
    end
    #back foward
    for i in row:-1:1
        x[i]=(a[i,col]-dot(a[i,i+1:row],x[i+1:row]))/a[i,i]
    end
    println("answer=",x)
    println("check ax-b=",ori_a*x-b)
end 

main(a,b,x)

Python版

import numpy as np

# Example1
a = np.matrix([[0, 4, 5, 2],
               [1, 0, 2, -6],
               [4, 1, 0, -2],
               [1, 7, 1, 0]], dtype='float64')
b = np.matrix([9, -3, 1, -3], dtype='float64').T
x = np.matrix([None]*4, dtype='float64').T

# for check sum
ori_a = a
# define extended coefficient matrix of eq ax=b
a = np.concatenate([a, b], axis=1)
(row, col) = a.shape
for j in range(row):
    # search pivot
    max_idx = j+np.argmax(abs(a[j:, j]))
    # swap!
    a[[j, max_idx]] = a[[max_idx, j]]
    # push forward
    for i in range(j+1, row):
        p = -a[i, j] / a[j, j]
        a[i] += p*a[j]
print('extended coefficient matrix applied gauss method \n a={}'.format(a))
# back forward
for i in range(row)[::-1]:
    x[i] = (a[i, -1]-np.dot(a[i, i+1:row], x[i+1:row]))/a[i, i]
print('x={}'.format(x))
# confirm
print("ax-b\n={}".format(ori_a@x-b))

すでに from numpy import * をしている状態でさぁコードを書いていきましょうという雰囲気で書けるんだなと感じました．型指定がめんどうですけれどね・・・．( 一一)