# islerles.jl -- numerical energy growth in a conservative system using Printf, LinearAlgebra logYs=zeros(10); loghs=zeros(10) function f(t,y) r=[y[2]*y[3]*sin(t)-y[1]*y[2]*y[3], -y[1]*y[3]*sin(t)+1/20*y[1]*y[3], y[1]^2*y[2]-1/20*y[1]*y[2]] return r end function euler(t,y,h) return y+h*f(t,y) end function rk2(t,y,h) k1=f(t,y) k2=f(t,y+h*2/3*k1) return y+h*(1/4*k1+3/4*k2) end function solve(t0,y0,h,n) yn=copy(y0) for j=1:n tn=t0+(j-1)*h yn=rk2(tn,yn,h) end return yn end function main() t0=0; T=1 y0=[1.0,-1,1] yold=zeros(3) @printf("%5s %7s %14s %14s %14s %14s\n", "p","n","h","E(1)","dE(1)","error") for p=6:16 n=2^p h=T/n yn=solve(t0,y0,h,n) @printf("%5d %7d %14.7e %14.7e %14.7e %14.7e\n", p,n,h,yn'*yn,yn'*yn-y0'*y0,norm(yold-yn)) if p>6 t=log(norm(yold-yn)) global logYs[p-6]=t global loghs[p-6]=log(2*h) end yold=copy(yn) end end main()