I'm doing a problem involving using numeric methods to solve a differential equation, and thought this might be a good opportunity to learn some basic Haskell. I have the following recurrence relation:
and the initial condition u(x, 0) = x^2
. I translated these into Haskell like so (putting in appropriate values for a,b,c,h and k from the specific problem, and noting that u_ij is defined as u(i*h, j*k)
):
u :: (Floating a, Eq a) => a -> a -> a
u x 0 = x*x
u x t = a*k / b*h * (u (x-h) (t-k))
+ (1 - (3*k/2*h))*(u x (t-k))
+ k/b * cos x
where
a = 3
b = 2
k = 0.1
h = 0.2
main = putStrLn (show (u 1 0.5))
This appears to run indefinitely. My best guess as to why is that floating point representation means the u x 0
pattern never actually matches. The way I'm used to dealing with this in other languages is to check if the absolute difference between the values is within some suitable epsilon, but this doesn't appear to be available for pattern matching. So it seems that floating point and pattern matching are fundamentally incompatible. Is this the likely problem, and if so, is there a canonical way to avoid one or the other in a situation like this?
It turns out the solution was to read a little further into the tutorial and use guards instead of pattern matching; they make checking against an epsilon quite trivial:
u :: (Floating a, Ord a) => a -> a -> a
u x t
| abs t <= 0.1 = x*x
| otherwise = a*k / b*h * (u (x-h) (t-k))
+ (1 - (3*k/2*h))*(u x (t-k))
+ k/b * cos x
where
a = 3
b = 2
k = 0.1
h = 0.2
main = putStrLn (show (u 1 0.5))
Gives an answer very quickly.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With