julia all the code of the paper, python still missing a notebook.

Notebooks/julia/DeC.jl

+#using Plots
+using FastGaussQuadrature, Compat
+using LinearAlgebra
+using Polynomials
+
+"""Different models"""
+
+function lagrange_basis(X,t,i)
+    idxs = eachindex(X)
+    if t in X
+        if t == X[i]
+            return 1
+        else
+            return 0
+        end
+    end
+    return prod((t-X[j])/(X[i]-X[j]) for j in idxs if j != i)
+end
+
+            
+function equispaced(order)
+    nodes=range(-1,stop=1,length=order)
+    w=zeros(order)
+    for k=1:order
+        zz=Polynomial([1.])
+        for s=1:order
+            if s == k
+                continue
+            else
+                zz=zz* Polynomial([-nodes[s],1])
+            end
+        end
+        zz=zz/zz(nodes[k])
+        intzz=integrate(zz)
+        w[k]=intzz(1)-intzz(-1)
+    end
+    return nodes, w
+end
+     
+
+function get_nodes(order,nodes_type)
+    if nodes_type=="equispaced"
+        nodes, w = equispaced(order)
+    elseif nodes_type=="gaussLobatto"
+        nodes, w = gausslobatto(order)
+    elseif nodes_type=="gaussLegendre"
+        nodes, w = gausslegendre(order)
+    end
+    w = map( x-> x*0.5, w)
+    nodes = nodes*0.5.+0.5    
+    return nodes, w
+end
+
+
+function compute_theta_DeC(order, nodes_type)
+    
+    nodes, w = get_nodes(order,nodes_type)
+    int_nodes, int_w = get_nodes(order,"gaussLobatto")
+    # generate theta coefficients 
+    theta = zeros(order,order)
+    beta_m =zeros(order)
+    for m in 1:order
+        beta_m[m] =nodes[m]
+        nodes_m = int_nodes*(nodes[m])
+        w_m = int_w*(nodes[m])
+        for r in 1:order
+            theta[r,m] = sum([lagrange_basis(nodes,nodes_m[k],r)*w_m[k] for k in 1:order])
+        end
+    end
+    return theta, beta_m
+    
+end
+            
+function compute_theta_DeC_subtimesteps(order, nodes_type)
+    
+    nodes, w = get_nodes(order,nodes_type)
+    int_nodes, int_w = get_nodes(order,"gaussLobatto")
+    # generate theta coefficients 
+    theta = zeros(order,order)
+    beta_m =zeros(order)
+    for m in 1:order
+        if m==1
+            beta_m[m] = nodes[1]
+            nodes_m=int_nodes*beta_m[m] 
+            w_m = int_w*beta_m[m]
+        else
+            beta_m[m] = (nodes[m]-nodes[m-1])
+            nodes_m = int_nodes*beta_m[m] .+nodes[m-1]
+            w_m = int_w*beta_m[m]
+        end
+        for r in 1:order
+            theta[r,m] = sum([lagrange_basis(nodes,nodes_m[k],r)*w_m[k] for k in 1:order])
+        end
+    end
+    return theta, beta_m
+    
+end
+
+
+""" Computing the matrix of the values T_k^m = int_{t^0}^{t^m} varphi_k"""
+function calculate_theta(M_sub::Int, distribution)
+  if distribution=="equispaced"
+    nodes=range(0,stop=1, length=M_sub+1)
+  elseif distribution=="gaussLobatto"
+    nodes,p = gausslobatto(M_sub+1)
+    nodes=nodes*0.5.+0.5
+  end
+  theta=compute_theta(nodes)
+  return theta
+end
+
+
+#Dec_correction-
+#Func is the function of the ode
+#N_step -> Number of subtimesteps
+#Corr-> Number of Corrections
+#y_0 Initial condition
+# t_0 Start value
+# t_end End value
+# Maybe logical variable for checking for the different solvers or the distributaiton
+# of subtimesteps up to some order
+#Distribution is set equidistand or different order
+
+""" Dec_correction constant time step"""
+function dec(func, tspan, y_0, M_sub::Int, K_corr::Int, distribution)
+
+   N_time=length(tspan)
+   dim=length(y_0)
+   U=zeros(dim, N_time)
+   u_p=zeros(dim, M_sub+1)
+   u_a=zeros(dim, M_sub+1)
+   u_help=zeros(dim)
+   func_corr=zeros(dim,M_sub+1)
+   Theta, beta = compute_theta_DeC(M_sub+1,distribution)
+   U[:,1]=y_0
+   for it=2: N_time
+        delta_t=(tspan[it]-tspan[it-1])
+        for m=1:M_sub+1
+         u_a[:,m]=U[:,it-1]
+         u_p[:,m]=U[:,it-1]
+        end
+        for k=2:K_corr+1
+            u_p=copy(u_a)
+            for r=1:M_sub+1
+                func_corr[:,r]=func(u_p[:,r])
+            end
+            for m=2:M_sub+1
+            u_a[:,m]= u_a[:,1]+delta_t*sum(Theta[r,m]*func_corr[:,r] for r in 1:M_sub+1)
+            end
+        end
+        U[:,it]=u_a[:,M_sub+1]
+    end
+    return tspan, U
+end
+            
+function dec_subtimesteps(func, tspan, y_0, M_sub::Int, K_corr::Int, distribution)
+
+   N_time=length(tspan)
+   dim=length(y_0)
+   U=zeros(dim, N_time)
+   u_p=zeros(dim, M_sub+1)
+   u_a=zeros(dim, M_sub+1)
+   u_help=zeros(dim)
+   func_corr=zeros(dim,M_sub+1)
+   Theta, beta = compute_theta_DeC_subtimesteps(M_sub+1,distribution)
+   U[:,1]=y_0
+   for it=2: N_time
+        delta_t=(tspan[it]-tspan[it-1])
+        for m=1:M_sub+1
+         u_a[:,m]=U[:,it-1]
+         u_p[:,m]=U[:,it-1]
+        end
+        for k=2:K_corr+1
+            u_p=copy(u_a)
+            for r=1:M_sub+1
+                func_corr[:,r]=func(u_p[:,r])
+            end
+            for m=2:M_sub+1
+            u_a[:,m]= u_a[:,m-1]+delta_t*sum(Theta[r,m]*func_corr[:,r] for r in 1:M_sub+1)
+            end
+        end
+        U[:,it]=u_a[:,M_sub+1]
+    end
+    return tspan, U
+end
+
+""" Dec_correction constant time step""" 
+function energy_relaxed_dec(func, t0, t_end, dt0, y_0, M_sub::Int, K_corr::Int, distribution)
+
+   dim=length(y_0)
+   U=zeros(dim,1)
+   U[:,1]=y_0
+   time=t0
+   tSpan=t0
+   gammas=1.0
+   it=1
+   u_p=zeros(dim, M_sub+1)
+   u_a=zeros(dim, M_sub+1)
+   u_help=zeros(dim)
+   func_corr=zeros(dim,M_sub+1)
+   Theta, beta = compute_theta_DeC(M_sub+1,distribution)
+   while (time <t_end)
+        it=it+1
+        delta_t=dt0
+        for m=1:M_sub+1
+         u_a[:,m]=U[:,it-1]
+         u_p[:,m]=U[:,it-1]
+        end
+        for k=2:K_corr+1
+            u_p=copy(u_a)
+            for r=1:M_sub+1
+                func_corr[:,r]=func(u_p[:,r])
+            end
+            for m=2:M_sub+1
+            u_a[:,m]= u_a[:,1]+delta_t*sum(Theta[r,m]*func_corr[:,r] for r in 1:M_sub+1)
+            end
+        end
+        udiff=u_a[:,M_sub+1]-U[:,it-1]
+        gamma=-2*(udiff'*U[:,it-1])/(udiff'*udiff)
+        delta_t_relaxed=gamma*delta_t
+        time= time+delta_t_relaxed
+        tSpan=[tSpan time]
+        gammas=[gammas gamma]
+        uRelaxed = U[:,it-1]+gamma*udiff
+        U=[U uRelaxed]
+    end
+    return tSpan, U, gammas
+end
+
+""" Dec_correction constant time step""" 
+function energy_relaxed_dec_dissipative(func, t0, t_end, dt0, y_0, M_sub::Int, K_corr::Int, distribution)
+
+   dim=length(y_0)
+   U=zeros(dim,1)
+   U[:,1]=y_0
+   time=t0
+   tSpan=t0
+   gammas=1.0
+   it=1
+   u_p=zeros(dim, M_sub+1)
+   u_a=zeros(dim, M_sub+1)
+   u_help=zeros(dim)
+   func_corr=zeros(dim,M_sub+1)
+   Theta, beta = compute_theta_DeC(M_sub+1,distribution)
+   while (time <t_end)
+        it=it+1
+        delta_t=dt0
+        for m=1:M_sub+1
+         u_a[:,m]=U[:,it-1]
+         u_p[:,m]=U[:,it-1]
+        end
+        for k=2:K_corr+1
+            u_p=copy(u_a)
+            for r=1:M_sub+1
+                func_corr[:,r]=func(u_p[:,r])
+            end
+            for m=2:M_sub+1
+                u_a[:,m]= u_a[:,1]+delta_t*sum(Theta[r,m]*func_corr[:,r] for r in 1:M_sub+1)
+            end
+        end
+        udiff=u_a[:,M_sub+1]-U[:,it-1]
+        gamma=2*delta_t*sum(Theta[r,M_sub+1]*func_corr[:,r]'*(u_p[:,r]-U[:,it-1])  for r in 1:M_sub+1)/(udiff'*udiff)
+        delta_t_relaxed=gamma*delta_t
+        time= time+delta_t_relaxed
+        tSpan=[tSpan time]
+        gammas=[gammas gamma]
+        uRelaxed = U[:,it-1]+gamma*udiff
+        U=[U uRelaxed]
+    end
+    return tSpan, U, gammas
+end
+      
+
+function dec_stiff(func, jac_stiff, tspan, y_0, M_sub::Int, K_corr::Int, distribution)
+   verbose=0
+   N_time=length(tspan)
+   dim=length(y_0)
+   U=zeros(dim, N_time)
+   u_p=zeros(dim, M_sub+1)
+   u_a=zeros(dim, M_sub+1)
+   u_help=zeros(dim)
+   func_corr=zeros(dim,M_sub+1)
+   Theta, beta = compute_theta_DeC(M_sub+1,distribution)
+   invJac=zeros(M_sub+1,dim,dim)
+   U[:,1]=y_0
+   for it=2: N_time
+        delta_t=(tspan[it]-tspan[it-1])
+        for m=1:M_sub+1
+         u_a[:,m]=U[:,it-1]
+         u_p[:,m]=U[:,it-1]
+        end
+        SS=jac_stiff(u_p[:,1])
+        for m=2:M_sub+1
+            invJac[m,:,:]=inv(Matrix(I,dim,dim).+ delta_t*beta[m] .*SS)
+        end
+        for k=2:K_corr+1
+            u_p=copy(u_a)
+            for r=1:M_sub+1
+                func_corr[:,r]=func(u_p[:,r])
+            end
+            for m=2:M_sub+1
+                u_a[:,m]= u_p[:,m]+delta_t.*invJac[m,:,:]*(-(u_p[:,m]-u_p[:,1])./delta_t+sum(Theta[r,m]*func_corr[:,r] for r in 1:M_sub+1))
+            end
+        if verbose==1
+            println("dt $(delta_t)")
+            println("up $(u_p)")
+            println("jac $(jac_stiff(u_p[:,1]))")
+            println("invJac $(invJac)")
+        end
+        end
+        U[:,it]=u_a[:,M_sub+1]
+    end
+    return tspan, U
+end
+
+            
+                        
+function dec_stiff_subtimesteps(func, jac_stiff, tspan, y_0, M_sub::Int, K_corr::Int, distribution)
+   verbose=0
+   N_time=length(tspan)
+   dim=length(y_0)
+   U=zeros(dim, N_time)
+   u_p=zeros(dim, M_sub+1)
+   u_a=zeros(dim, M_sub+1)
+   u_help=zeros(dim)
+   func_corr=zeros(dim,M_sub+1)
+   Theta, beta = compute_theta_DeC_subtimesteps(M_sub+1,distribution)
+   S = zeros(dim,dim)
+   invJac=zeros(M_sub+1,dim,dim)
+   U[:,1]=y_0
+   for it=2: N_time
+        delta_t=(tspan[it]-tspan[it-1])
+        for m=1:M_sub+1
+         u_a[:,m]=U[:,it-1]
+         u_p[:,m]=U[:,it-1]
+        end
+        SS=jac_stiff(u_p[:,1])
+        for m=1:M_sub+1
+            invJac[m,:,:]=inv(Matrix(I,dim,dim) .+ delta_t*beta[m] .*SS)
+        end
+        for k=2:K_corr+1
+            u_p=copy(u_a)
+            for r=1:M_sub+1
+                func_corr[:,r]=func(u_p[:,r])
+            end
+            for m=2:M_sub+1
+                L2=(u_p[:,m]-u_a[:,m-1])-delta_t.*sum(Theta[r,m]*func_corr[:,r] for r in 1:M_sub+1) #
+                u_a[:,m]= u_p[:,m]+ invJac[m,:,:]*( - L2) #u_a[:,m-1]-u_p[:,m-1]
+                #u_a[:,m]= invJac[m,:,:]*(u_a[:,m-1] +delta_t*beta[m] .*SS*u_p[:,m] -L2)
+            end
+        if verbose==1
+                            println("beta $(beta)")
+                            println("theta $(Theta)")
+            println("dt $(delta_t)")
+            println("up $(u_p)")
+            println("jac $(jac_stiff(u_p[:,1]))")
+            println("invJac $(invJac)")
+        end
+        end
+        U[:,it]=u_a[:,M_sub+1]
+    end
+    return tspan, U
+end
+
+

Notebooks/julia/Lagrange_polynomials.jl

+using SymPy
+
+import SymPy
+
+
+function generate_poly(s, nodes,deg)
+#Generates Lagrangian polynomial in symbolic
+
+aux=[]
+phi=[]
+
+for k=1:deg+1
+    append!(aux, 1);
+    for j=1:deg+1
+        if j != k
+            aux[k]=aux[k].*(s-nodes[j]);
+        end
+    end
+end
+for j=1:deg+1
+    aux[j]= aux[j]/SymPy.N(subs(aux[j],s,nodes[j]));
+end
+return aux
+end
+
+function compute_theta(nodes)
+# Compute the integral of lagrange polynomials from 0 to interpolation point theta[r,m]= \int_{t^0}^{t^m} \phi_r
+s=Sym("s");
+deg=length(nodes)-1;
+poly=generate_poly(s,nodes,deg)
+
+theta=zeros(deg+1, deg+1)
+
+for r=1:deg+1
+    for m=1:deg+1
+        theta[r,m]=SymPy.integrate(poly[r],(s,0,nodes[m]))
+    end
+end
+
+return theta
+
+end

Notebooks/python/DeC.py

+import numpy as np
+import matplotlib.pyplot as plt
+from numpy.polynomial.legendre import leggauss
+from quadr import lglnodes,equispaced
+def lagrange_basis(nodes,x,k):
+    y=np.zeros(x.size)
+    for ix, xi in enumerate(x):
+        tmp=[(xi-nodes[j])/(nodes[k]-nodes[j])  for j in range(len(nodes)) if j!=k]
+        y[ix]=np.prod(tmp)
+    return y
+
+def get_nodes(order,nodes_type):
+    if nodes_type=="equispaced":
+        nodes,w = equispaced(order)
+    elif nodes_type == "gaussLegendre":
+        nodes,w = leggauss(order)
+    elif nodes_type == "gaussLobatto":
+        nodes, w = lglnodes(order-1,10**-15)
+    nodes=nodes*0.5+0.5
+    w = w*0.5
+    return nodes, w
+        
+def compute_theta_DeC(order, nodes_type):
+    nodes, w = get_nodes(order,nodes_type)
+    int_nodes, int_w = get_nodes(order,"gaussLobatto")
+    # generate theta coefficients 
+    theta = np.zeros((order,order))
+    beta = np.zeros(order)
+    for m in range(order):
+        beta[m] = nodes[m]
+        nodes_m = int_nodes*(nodes[m])
+        w_m = int_w*(nodes[m])
+        for r in range(order):
+            theta[r,m] = sum(lagrange_basis(nodes,nodes_m,r)*w_m)
+    return theta, beta
+
+
+def compute_RK_from_DeC(M_sub,K_corr,nodes_type):
+    order=M_sub+1;
+    [theta,beta]=compute_theta_DeC(order,nodes_type)
+    bar_beta=beta[1:]  # M_sub
+    bar_theta=theta[:,1:].transpose() # M_sub x (M_sub +1)
+    theta0= bar_theta[:,0]  # M_sub x 1
+    bar_theta= bar_theta[:,1:] #M_sub x M_sub
+    A=np.zeros((M_sub*(K_corr-1)+1,M_sub*(K_corr-1)+1))  # (M_sub x K_corr +1)^2
+    b=np.zeros(M_sub*(K_corr-1)+1)
+    c=np.zeros(M_sub*(K_corr-1)+1)
+
+    c[1:M_sub+1]=bar_beta
+    A[1:M_sub+1,0]=bar_beta
+    for k in range(1,K_corr-1):
+        r0=1+M_sub*k
+        r1=1+M_sub*(k+1)
+        c0=1+M_sub*(k-1)
+        c1=1+M_sub*(k)
+        c[r0:r1]=bar_beta
+        A[r0:r1,0]=theta0
+        A[r0:r1,c0:c1]=bar_theta
+    b[0]=theta0[-1]
+    b[-M_sub:]=bar_theta[M_sub-1,:]
+    return A,b,c
+
+
+def dec(func, tspan, y_0, M_sub, K_corr, distribution):
+    N_time=len(tspan)
+    dim=len(y_0)
+    U=np.zeros((dim, N_time))
+    u_p=np.zeros((dim, M_sub+1))
+    u_a=np.zeros((dim, M_sub+1))
+    rhs= np.zeros((dim,M_sub+1))
+    Theta, beta = compute_theta_DeC(M_sub+1,distribution)
+    U[:,0]=y_0
+    for it in range(1, N_time):
+        delta_t=(tspan[it]-tspan[it-1])
+        for m in range(M_sub+1):
+            u_a[:,m]=U[:,it-1]
+            u_p[:,m]=U[:,it-1]
+        for k in range(1,K_corr+1):
+            u_p=np.copy(u_a)
+            for r in range(M_sub+1):
+                rhs[:,r]=func(u_p[:,r])
+            for m in range(1,M_sub+1):
+                u_a[:,m]= U[:,it-1]+delta_t*sum([Theta[r,m]*rhs[:,r] for r in range(M_sub+1)])
+        U[:,it]=u_a[:,M_sub]
+    return tspan, U
+            
+def decImplicit(func,jac_stiff, tspan, y_0, M_sub, K_corr, distribution):
+    N_time=len(tspan)
+    dim=len(y_0)
+    U=np.zeros((dim, N_time))
+    u_p=np.zeros((dim, M_sub+1))
+    u_a=np.zeros((dim, M_sub+1))
+    u_help= np.zeros(dim)
+    rhs= np.zeros((dim,M_sub+1))
+    Theta, beta = compute_theta_DeC(M_sub+1,distribution)
+    invJac=np.zeros((M_sub+1,dim,dim))
+    U[:,0]=y_0
+    for it in range(1, N_time):
+        delta_t=(tspan[it]-tspan[it-1])
+        for m in range(M_sub+1):
+            u_a[:,m]=U[:,it-1]
+            u_p[:,m]=U[:,it-1]
+        SS=jac_stiff(u_p[:,0])
+        for m in range(1,M_sub+1):
+            invJac[m,:,:]=np.linalg.inv(np.eye(dim) - delta_t*beta[m]*SS)
+        for k in range(1,K_corr+1):
+            u_p=np.copy(u_a)
+            for r in range(M_sub+1):
+                rhs[:,r]=func(u_p[:,r])
+            for m in range(1,M_sub+1):
+                u_a[:,m]= u_p[:,m]+delta_t*np.matmul(invJac[m,:,:],\
+                (-(u_p[:,m]-u_p[:,0])/delta_t\
+                 +sum([Theta[r,m]*rhs[:,r] for r in range(M_sub+1)])))
+        U[:,it]=u_a[:,M_sub]
+    return tspan, U
+
+
+
+def decMPatankar(prod_dest, rhs, tspan, y_0, M_sub, K_corr, distribution):
+    N_time=len(tspan)
+    dim=len(y_0)
+    U=np.zeros((dim, N_time))
+    u_p=np.zeros((dim, M_sub+1))
+    u_a=np.zeros((dim, M_sub+1))
+    prod_p = np.zeros((dim,dim,M_sub+1))
+    dest_p = np.zeros((dim,dim,M_sub+1))
+    rhs_p= np.zeros((dim,M_sub+1))
+    Theta, beta = compute_theta_DeC(M_sub+1,distribution)
+    U[:,0]=y_0
+    for it in range(1, N_time):
+        delta_t=(tspan[it]-tspan[it-1])
+        for m in range(M_sub+1):
+            u_a[:,m]=U[:,it-1]
+            u_p[:,m]=U[:,it-1]
+        for k in range(1,K_corr+1):
+            u_p=np.copy(u_a)
+            for r in range(M_sub+1):
+                prod_p[:,:,r], dest_p[:,:,r]=prod_dest(u_p[:,r])
+                rhs_p[:,r]=rhs(u_p[:,r])
+            for m in range(1,M_sub+1):
+                u_a[:,m]= patankar_type_dec(prod_p,dest_p,rhs_p,delta_t,m,M_sub,Theta,u_p,dim)
+        U[:,it]=u_a[:,M_sub]
+    return tspan, U
+
+
+def patankar_type_dec(prod_p,dest_p,rhs_p,delta_t,m,M_sub,Theta,u_p,dim):
+    mass= np.eye(dim)
+    RHS= u_p[:,0]
+    for i in range(dim):
+        for r in range(M_sub+1):
+            RHS[i]=RHS[i]+delta_t*Theta[r,m]*rhs_p[i,r]
+            if Theta[r,m]>0:
+                for j in range(dim):
+                    mass[i,j]=mass[i,j]-delta_t*Theta[r,m]*(prod_p[i,j,r]/u_p[j,m])
+                    mass[i,i]=mass[i,i]+ delta_t*Theta[r,m]*(dest_p[i,j,r]/u_p[i,m])
+            elif Theta[r,m]<0:
+                for j in range(dim):
+                    mass[i,i]=mass[i,i]- delta_t*Theta[r,m]*(prod_p[i,j,r]/u_p[i,m])
+                    mass[i,j]=mass[i,j]+ delta_t*Theta[r,m]*(dest_p[i,j,r]/u_p[j,m])
+    return np.linalg.solve(mass,RHS)

Notebooks/python/ODEproblems.py

+import numpy as np
+
+## Linear scalar Dahlquist's equation
+def linear_scalar_flux(u,t=0,k_coef=10):
+    ff=np.zeros(np.shape(u))
+    ff[0]= -k_coef*u[0]
+    return ff
+
+def linear_scalar_exact_solution(u0,t,k_coef=10):
+    return np.array([np.exp(-k_coef*u0[0]*t)])