28 lines
952 B
Ruby
28 lines
952 B
Ruby
def calc_rk4(f)
|
|
return ->(t,y,dt){
|
|
->(dy1 ){
|
|
->(dy2 ){
|
|
->(dy3 ){
|
|
->(dy4 ){ ( dy1 + 2*dy2 + 2*dy3 + dy4 ) / 6 }.call(
|
|
dt * f.call( t + dt , y + dy3 ))}.call(
|
|
dt * f.call( t + dt/2, y + dy2/2 ))}.call(
|
|
dt * f.call( t + dt/2, y + dy1/2 ))}.call(
|
|
dt * f.call( t , y ))}
|
|
end
|
|
|
|
TIME_MAXIMUM, WHOLE_TOLERANCE = 10.0, 1.0e-5
|
|
T_START, Y_START, DT = 0.0, 1.0, 0.10
|
|
|
|
def my_diff_eqn(t,y) ; t * Math.sqrt(y) ; end
|
|
def my_solution(t ) ; (t**2 + 4)**2 / 16 ; end
|
|
def find_error(t,y) ; (y - my_solution(t)).abs ; end
|
|
def is_whole?(t ) ; (t.round - t).abs < WHOLE_TOLERANCE ; end
|
|
|
|
dy = calc_rk4( ->(t,y){my_diff_eqn(t,y)} )
|
|
|
|
t, y = T_START, Y_START
|
|
while t <= TIME_MAXIMUM
|
|
printf("y(%4.1f)\t= %12.6f \t error: %12.6e\n",t,y,find_error(t,y)) if is_whole?(t)
|
|
t, y = t + DT, y + dy.call(t,y,DT)
|
|
end
|