Custom linear system solver

New to Julia
1

Hi there,

I am new to julia. I am trying to write a function to solve tridiagonal matrix system by Thomas Algorithm

Here is my code

function TDMA(a,b,c,D)
#  Thomas Algorithm
#  Vectors a,b,c have the same length
#  Should have a[1]=0 and c[n]=0
#  A*x = D
#  [b[1]   c[1]               ] [ x[1] ]   [ D[1] ]
#  [a[2]   b[2]   c[2]        ] [ x[2] ]   [ D[2] ]
#  [    ...   ...   ...       ] [ ...  ] = [ ...  ]
#  [      a[n-1] b[n-1] c[n-1]] [x[n-1]]   [D[n-1]]
#  [              a[n]   b[n] ] [ x[n] ]   [ D[n] ]  
x = zeros(length(D),1)
d = D
n = length(D)
w = zeros(n,1)
    for i = 2:n
        w[i] = a[i]/b[i-1]
        b[i] = b[i] - w[i]*c[i-1]
        d[i] = d[i] - w[i]*d[i-1]
    end
    x[n] = d[n]/b[n]
    for k = 1:n-1
        x[n-k] = (d[n-k]-c[n-k]*x[n-k+1])/b[n-k]
    end
    x
end

I am testing it by this example

A = 
[1 2 0;
 4 5 6;
 0 1 2]
Y = [1; 3; 4]

with

a = [0 4 1]
b = [1 5 2]
c = [2 6 0]

to be the arguments.

The build-in function

A\Y

gives the solution

 -3.333333333333333 
  2.1666666666666665
  0.9166666666666666

Howevery, when I try to use my TDMA to solve the equation, it comes error

TDMA(a,b,c,Y)

InexactError: Int64(4.545454545454545)

Stacktrace:
 [1] Type at ./float.jl:703 [inlined]
 [2] convert at ./number.jl:7 [inlined]
 [3] setindex! at ./array.jl:767 [inlined]
 [4] TDMA(::Array{Int64,2}, ::Array{Int64,2}, ::Array{Int64,2}, ::Array{Int64,1}) at ./In[13]:17
 [5] top-level scope at In[14]:1

Can anyone help me with this?

Thanks

1 Like
2

Hi

The problem is the assignment of a floating point number to a vector of integers in the lines b[i] = b[i] - w[i]*c[i-1] and d[i] = d[i] - w[i]*d[i-1].

You see a similar error when you try to do

julia> b = [1, 2, 3]
3-element Array{Int64,1}:
 1
 2
 3

julia> b[1] = .3
ERROR: InexactError: Int64(0.3)
Stacktrace:
 [1] Type at ./float.jl:703 [inlined]
 [2] convert at ./number.jl:7 [inlined]
 [3] setindex!(::Array{Int64,1}, ::Float64, ::Int64) at ./array.jl:767
 [4] top-level scope at none:0

The solution is to either call your function with arrays of floats, e.g. by writing b = [1., 5., 2.], or convert these vectors inside you function with e.g. b = float.(b) and d = float.(D).
Note also that it is common to append ! to function names that change the arguments.

2 Likes
3

It works now. Thank you!

4

Note that there is a Tridiagonal type in LinearAlgebra by the way: Linear Algebra · The Julia Language. \ is probably overloaded for this type.

2 Likes
5

Yes but to me, it is slower than thomas. However, I think TDMA is numerically less stable

2 Likes