-
-
Save stillyslalom/8f9e4ab1a2c86a6cb0d6b34ade18d5e1 to your computer and use it in GitHub Desktop.
IRK
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using Printf, LinearAlgebra, StaticArrays | |
s3o6=sqrt(3)/6 | |
const a= SA[1/4 1/4-s3o6; 1/4+s3o6 1/4] | |
const b= SA[1/2,1/2] | |
const c= SA[1/2-s3o6,1/2+s3o6] | |
function euler(t,y,h) | |
return y+h*fF(t,y) | |
end | |
function rk2(t,y,h) | |
k1=fF(t,y) | |
k2=fF(t,y+h*2/3*k1) | |
return y+h*(1/4*k1+3/4*k2) | |
end | |
@inbounds function fF(t,y) | |
SA[ | |
sin(t)*y[2]*y[3] - y[1]*y[2]*y[3], | |
(y[1]*y[3])/20 - sin(t)*y[1]*y[3], | |
-(y[1]*y[2])/20 + y[1]^2*y[2] | |
] | |
end | |
@inbounds function fG(t,y,h,k) | |
SA[ | |
k[1] - sin(t + h*c[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]) + (h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]), | |
k[2] - sin(t + h*c[1])*(-(h*(a[1, 1]*k[1] + a[1, 2]*k[4])) - y[1])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]) - ((h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]))/20, | |
k[3] + ((h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2]))/20 - (h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])^2*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2]), | |
k[4] - sin(t + h*c[2])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]) + (h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]), | |
k[5] - sin(t + h*c[2])*(-(h*(a[2, 1]*k[1] + a[2, 2]*k[4])) - y[1])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]) - ((h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]))/20, | |
k[6] + ((h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2]))/20 - (h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])^2*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2]), | |
] | |
end | |
@inbounds function fDG(t,y,h,k) | |
DG= @MMatrix zeros(6,6) | |
DG[1, 1] = 1 + h*a[1, 1]*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]) | |
DG[1, 2] = -(sin(t + h*c[1])*h*a[1, 1]*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3])) + h*a[1, 1]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]) | |
DG[1, 3] = -(sin(t + h*c[1])*h*a[1, 1]*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])) + h*a[1, 1]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2]) | |
DG[1, 4] = h*a[1, 2]*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]) | |
DG[1, 5] = -(sin(t + h*c[1])*h*a[1, 2]*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3])) + h*a[1, 2]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]) | |
DG[1, 6] = -(sin(t + h*c[1])*h*a[1, 2]*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])) + h*a[1, 2]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2]) | |
DG[2, 1] = -(h*a[1, 1]*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]))/20 + sin(t + h*c[1])*h*a[1, 1]*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]) | |
DG[2, 2] = 1 | |
DG[2, 3] = -(sin(t + h*c[1])*h*a[1, 1]*(-(h*(a[1, 1]*k[1] + a[1, 2]*k[4])) - y[1])) - (h*a[1, 1]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1]))/20 | |
DG[2, 4] = -(h*a[1, 2]*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]))/20 + sin(t + h*c[1])*h*a[1, 2]*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]) | |
DG[2, 5] = 0 | |
DG[2, 6] = -(sin(t + h*c[1])*h*a[1, 2]*(-(h*(a[1, 1]*k[1] + a[1, 2]*k[4])) - y[1])) - (h*a[1, 2]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1]))/20 | |
DG[3, 1] = (h*a[1, 1]*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2]))/20 - 2*h*a[1, 1]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2]) | |
DG[3, 2] = (h*a[1, 1]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1]))/20 - h*a[1, 1]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])^2 | |
DG[3, 3] = 1 | |
DG[3, 4] = (h*a[1, 2]*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2]))/20 - 2*h*a[1, 2]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2]) | |
DG[3, 5] = (h*a[1, 2]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1]))/20 - h*a[1, 2]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])^2 | |
DG[3, 6] = 0 | |
DG[4, 1] = h*a[2, 1]*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]) | |
DG[4, 2] = -(sin(t + h*c[2])*h*a[2, 1]*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3])) + h*a[2, 1]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]) | |
DG[4, 3] = -(sin(t + h*c[2])*h*a[2, 1]*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])) + h*a[2, 1]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2]) | |
DG[4, 4] = 1 + h*a[2, 2]*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]) | |
DG[4, 5] = -(sin(t + h*c[2])*h*a[2, 2]*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3])) + h*a[2, 2]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]) | |
DG[4, 6] = -(sin(t + h*c[2])*h*a[2, 2]*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])) + h*a[2, 2]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2]) | |
DG[5, 1] = -(h*a[2, 1]*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]))/20 + sin(t + h*c[2])*h*a[2, 1]*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]) | |
DG[5, 2] = 0 | |
DG[5, 3] = -(sin(t + h*c[2])*h*a[2, 1]*(-(h*(a[2, 1]*k[1] + a[2, 2]*k[4])) - y[1])) - (h*a[2, 1]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1]))/20 | |
DG[5, 4] = -(h*a[2, 2]*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]))/20 + sin(t + h*c[2])*h*a[2, 2]*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]) | |
DG[5, 5] = 1 | |
DG[5, 6] = -(sin(t + h*c[2])*h*a[2, 2]*(-(h*(a[2, 1]*k[1] + a[2, 2]*k[4])) - y[1])) - (h*a[2, 2]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1]))/20 | |
DG[6, 1] = (h*a[2, 1]*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2]))/20 - 2*h*a[2, 1]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2]) | |
DG[6, 2] = (h*a[2, 1]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1]))/20 - h*a[2, 1]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])^2 | |
DG[6, 3] = 0 | |
DG[6, 4] = (h*a[2, 2]*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2]))/20 - 2*h*a[2, 2]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2]) | |
DG[6, 5] = (h*a[2, 2]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1]))/20 - h*a[2, 2]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])^2 | |
DG[6, 6] = 1 | |
return DG | |
end | |
@inbounds function irk(t,y,h) | |
k1=fF(t+h*c[1],rk2(t,y,h*c[1])) | |
k2=fF(t+h*c[2],rk2(t,y,h*c[2])) | |
Ki=[k1;k2] | |
lastone=false | |
for i=1:10 | |
Gi=fG(t,y,h,Ki) | |
DGi=fDG(t,y,h,Ki) | |
DK=-DGi\Gi | |
Ki+=DK | |
if lastone | |
break | |
end | |
if norm(DK)<=norm(Ki)*10e-10 | |
lastone=true | |
end | |
end | |
k1=Ki[1:3] | |
k2=Ki[4:6] | |
return y+h*(b[1]*k1+b[2]*k2) | |
end | |
function solve(t0,y0,h,n) | |
yn=copy(y0) | |
for j=1:n | |
tn=t0+(j-1)*h | |
yn=irk(tn,yn,h) | |
end | |
return yn | |
end | |
function main(args) | |
@printf("conjulia -- Energy Conserving IRK Scheme\n\n") | |
t0=0; T=100 | |
y0=[1.0,1,1] | |
if length(args)==3 | |
for j=1:3 | |
y0[j]=parse(Float64,args[j]) | |
end | |
end | |
println("y($t0)=$y0\n") | |
einit=y0'*y0 | |
yold=zeros(3) | |
@printf("%5s %7s %14s %14s %14s %14s\n", | |
"p","n","h","E(100)","dE(100)","error") | |
yn=copy(y0) | |
n=64 | |
for p=6:16 | |
h=T/n | |
yn=solve(t0,y0,h,n) | |
if p>6 | |
@printf("%5d %7d %14.7e %14.7e %14.7e %14.7e\n", | |
p,n,h,yn'*yn,yn'*yn-einit,norm(yold-yn)) | |
end | |
flush(stdout) | |
yold=copy(yn) | |
n*=2 | |
end | |
println("\ny($T)=$yn") | |
end | |
main(ARGS) | |
tsec=@elapsed main(ARGS) | |
println("\nTotal execution time ",tsec," seconds.") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment