<jupyter_start><jupyter_text>Two PDEs (diffusion/heat equation and an advection-decay equation)<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from nodepy import rk
import cvxpy as cp

import matplotlib.lines as mlines

import plot_fkt
plot_fkt.setup_plt()


import scipy.linalg as linalg
import timeit


#Diagonally Implicit methods:
BE = rk.loadRKM('BE').__num__()
SDIRK23 = rk.loadRKM('SDIRK23').__num__()
SDIRK34 = rk.loadRKM('SDIRK34').__num__()
SDIRK54 = rk.loadRKM('SDIRK54').__num__()
TR_BDF2 = rk.loadRKM('TR-BDF2').__num__()


be = rk.loadRKM('BE').__num__()


#Extrapolation method
ex2 = rk.extrap(2,'implicit euler').__num__()
ex3 = rk.extrap(3,'implicit euler').__num__()
ex4 = rk.extrap(4,'implicit euler').__num__()
ex5 = rk.extrap(5,'implicit euler').__num__()
ex6 = rk.extrap(6,'implicit euler').__num__()
ex8 = rk.extrap(8,'implicit euler').__num__()

bs5 = rk.loadRKM('BS5').__num__()
dp5 = rk.loadRKM('DP5').__num__()

from OrderCondition import *
from RKimple import *
import utils 

from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, mark_inset
from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
def draw_order(x0=0,x1=1,y0=0,p=1,ftext=1.2,text='above'):
    K= y0/(x0**p)
    plt.plot([x0,x1],[K*x0**p,K*x1**p],'ko-')
    
    xtext = 10**(0.5*(np.log10(x0)+np.log10(x1)))
    if text=='above':
        ytext = ftext*K*xtext**p
        plt.text(xtext,ytext,'$p='+str(p)+'$',horizontalalignment='right',verticalalignment='bottom')
    elif text=='below':
        ytext = (1/ftext)*K*xtext**p
        plt.text(xtext,ytext,'$p='+str(p)+'$',horizontalalignment='left',verticalalignment='top')
    else:
        print("text = 'above' or text = 'below'")
        
def show_legend_marker(order=[4],loc='upper center'):
    original = mlines.Line2D([], [], color='black', marker=plot_fkt.marker(-1), label='$b$ unaltered')
    handles = [original]
    for o in order:
        handles.append(mlines.Line2D([], [], color='black', marker=plot_fkt.marker(o),label='$b$ adapted, order $\\geq'+ str(o)+'$'))
    legend_marker = plt.legend(handles=handles,loc=loc)
    plt.gca().add_artist(legend_marker)<jupyter_output><empty_output><jupyter_text>Diffusion/heat equation<jupyter_code>ex3 = rk.extrap(3,'implicit euler').__num__()
b_hat_1 = [np.array([0,0,0,1/3,1/3,1/3])]
b_hat_orig = [ex3.b]
ex3.b_hat = {3:b_hat_orig,2:b_hat_orig,1:b_hat_1}
N=100
x = np.linspace(0,1,N)
dx = x[1]-x[0]
u0 = np.zeros_like(x)
u0[int(N/2)] = 1
u0[int(N/2-1)] = 1
#dt = 0.7*dx**2
dt = 0.001

A_heat = 1/dx**2 * (-2*np.diag(np.ones(N))+np.diag(np.ones(N-1),-1)+np.diag(np.ones(N-1),1))
solver = Solver(rkm = ex3,
               dt = dt,
               t_final = 0.01,
               b_fixed=False,
               tol_neg=1e-8,
               tol_change = np.inf,
               p = [3,2,1],
               theta = [1],
               solver = cp.MOSEK,
               solver_eqs=solver_Matrix,
               convex = False,
               LP_opts = {'reduce':True})

problem_ADR = Problem(f=A_heat,
                 u0=u0,
                 minval=0,
                 maxval=np.inf)

status_ref,t,u,b,KK = RK_integrate(solver=solver,problem=problem_ADR,verbose=True,dumpK=True)

t_ref = np.array(t)
u_ref = np.array(u).T
b_ref = np.array(b).T
utils.show_status(status_ref)
solver = Solver(rkm = ex3,
               dt = dt,
               t_final = 0.01,
               b_fixed=False,
               tol_neg=1e-8,
               tol_change = np.inf,
               p = [1],
               theta = [1],
               solver = cp.MOSEK,
               solver_eqs=solver_Matrix,
               convex=True,
               LP_opts = {'reduce':True})

status_con,t,u,b,KK = RK_integrate(solver=solver,problem=problem_ADR,verbose=True,dumpK=True)

t_con = np.array(t)
u_con = np.array(u).T
b_con = np.array(b).T
utils.show_status(status_con)
index = [0,1,5,-1]
x=np.linspace(-0.5,0.5,100)
fig = plt.figure(figsize=[6.4, 4])
ax = fig.subplots()
for i in index:
    ax.plot(x,u_ref[:,i],label='t='+str(round(t_ref[i],3)))
    
plt.grid()
plt.xlabel('$x$')
plt.ylabel('$u$')
plt.legend(loc=2)
#plt.title('Diffusion Free Adaptation $\Delta t ='+str(dt)+'$')

axins = zoomed_inset_axes(ax, 2.2, loc=1)  # zoom = 6
for i in index:
    axins.plot(x,u_ref[:,i])

# sub region of the original image
x1, x2, y1, y2 = -0.1, 0.1, -0.01, 0.3
axins.set_xlim(x1, x2)
axins.set_ylim(y1, y2)
# fix the number of ticks on the inset axes
axins.yaxis.get_major_locator().set_params(nbins=7)
axins.xaxis.get_major_locator().set_params(nbins=7)
axins.grid()

plt.setp(axins.get_xticklabels(), visible=False)
plt.setp(axins.get_yticklabels(), visible=False)

# draw a bbox of the region of the inset axes in the parent axes and
# connecting lines between the bbox and the inset axes area
mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5")
plt.savefig('Diff_Direct.pdf',bbox_inches = "tight")
index = [0,1,5,-1]
x=np.linspace(-0.5,0.5,100)
fig = plt.figure(figsize=[6.4, 4])
ax = fig.subplots()
for i in index:
    ax.plot(x,u_con[:,i],label='$t='+str(round(t_ref[i],3))+'$')
    
plt.grid()
plt.xlabel('$x$')
plt.ylabel('$u$')
plt.legend(loc=2)
#plt.title('Diffusion Convex Adaptation $\Delta t ='+str(dt)+'$')

axins = zoomed_inset_axes(ax, 2.2, loc=1)  # zoom = 6
for i in index:
    axins.plot(x,u_con[:,i])

# sub region of the original image
x1, x2, y1, y2 = -0.1, 0.1, -0.01, 0.3
axins.set_xlim(x1, x2)
axins.set_ylim(y1, y2)
# fix the number of ticks on the inset axes
axins.yaxis.get_major_locator().set_params(nbins=7)
axins.xaxis.get_major_locator().set_params(nbins=7)
axins.grid()

plt.setp(axins.get_xticklabels(), visible=False)
plt.setp(axins.get_yticklabels(), visible=False)

# draw a bbox of the region of the inset axes in the parent axes and
# connecting lines between the bbox and the inset axes area
mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5")
plt.savefig('Diff_Convex.pdf',bbox_inches = "tight")
b_org = ex3.b.copy()
b_org.shape=(len(b_org),1)
t = t_ref

order = np.array(status_ref['order'])

plt.figure(figsize=[6.4, 2])
norm = np.linalg.norm((b_ref-b_org),axis=0,ord=1)
plt.plot(t[1:],norm[1:],'-',color='C1')
mask=(norm[1:]==0)
plot_fkt.plot_markers(t[1:],norm[1:],mask,order[1:],'C1')
plt.ylim([-0.15,3.7])
plt.xlabel('$t$')
plt.ylabel(r'$\| \tilde{b} - b \|_1 $')
plt.grid()
plt.legend()
plt.savefig('b_Diff_Direct.pdf',bbox_inches = "tight")
b_org = ex3.b.copy()
b_org.shape=(len(b_org),1)
t = t_con

order = np.array(status_con['order'])

plt.figure(figsize=[6.4, 2])
norm = np.linalg.norm((b_con-b_org),axis=0,ord=1)
plt.plot(t[1:],norm[1:],'-',color='C1')
mask=(norm[1:]==0)
plot_fkt.plot_markers(t[1:],norm[1:],mask,order[1:],'C1')
plt.ylim([-0.15,3.7])
plt.xlabel('$t$')
plt.ylabel(r'$\| \tilde{b} - b \|_1 $')
plt.grid()
plt.legend()
plt.savefig('b_Diff_Convex.pdf',bbox_inches = "tight")<jupyter_output><empty_output><jupyter_text>Convergence<jupyter_code>dts = np.logspace(-5,-2,num=20)
dts[-2] = 0.005

t_end = 0.01

reference = np.zeros([len(u0),len(dts)])
for i in range(len(dts)):
    reference[:,i]= linalg.expm(t_end*A_heat)@u0
    

problem_Heat = Problem(f=A_heat,
                 u0=u0,
                 minval=0,
                 maxval=np.inf)

solver_be = Solver(rkm = be,
               dt = dt,
               t_final = t_end,
               b_fixed=True,
               tol_neg=1e-8,
               tol_change = 5,
               p = [3,2,1],
               theta = [1],
               solver = cp.MOSEK,
               solver_eqs=solver_Matrix,
               LP_opts = {'reduce':True})

solver_ex3_fix = Solver(rkm = ex3,
               dt = dt,
               t_final = t_end,
               b_fixed=True,
               tol_neg=1e-8,
               tol_change = 5,
               p = [3,2,1],
               theta = [1],
               solver = cp.MOSEK,
               solver_eqs=solver_Matrix,
               LP_opts = {'reduce':True})

solver_ex3_adp = Solver(rkm = ex3,
               dt = dt,
               t_final = t_end,
               b_fixed=False,
               tol_neg=1e-8,
               tol_change = 5,
               p = [3,2,1],
               theta = [1],
               solver = cp.MOSEK,
               solver_eqs=solver_Matrix,
               convex=False,
               LP_opts = {'reduce':True,'verbose':False})


solver_ex3_cvx = Solver(rkm = ex3,
               dt = dt,
               t_final = t_end,
               b_fixed=False,
               tol_neg=1e-8,
               tol_change = 5,
               p = [1],
               theta = [1],
               solver = cp.MOSEK,
               solver_eqs=solver_Matrix,
               convex=True,
               LP_opts = {'reduce':True,'verbose':False})

print('BE')
sol_be,err_be,change_be = utils.plot_convergence(problem_Heat,solver_be,dts,reference,error='rel',step = -1,
                              Params=dict())

print('ex3')
sol_ex3,err_ex3,change_ex3 = utils.plot_convergence(problem_Heat,solver_ex3_fix,dts,reference,error='rel',step = -1,
                              Params=dict())

print('ex3_convex')
sol_ex3c,err_ex3c,change_ex3c,order_ex3c = utils.plot_convergence(problem_Heat,solver_ex3_cvx,dts,reference,error='rel',step = -1,
                              Params=dict(),get_order=True)

print('ex3_adapted')
sol_ex3a,err_ex3a,change_ex3a,order_ex3a = utils.plot_convergence(problem_Heat,solver_ex3_adp,dts,reference,error='rel',step = -1,
                              Params=dict(),get_order=True)
np.savez('convergence_heat.npz', 
         dts=dts,
        sol_be=sol_be, err_be=err_be, change_be=change_be,
        sol_ex3 = sol_ex3, err_ex3 = err_ex3, change_ex3 = change_ex3,
        sol_ex3c=sol_ex3c, err_ex3c=err_ex3c, change_ex3c=change_ex3c, order_ex3c=order_ex3c,
        sol_ex3a=sol_ex3a, err_ex3a=err_ex3a, change_ex3a=change_ex3a, order_ex3a=order_ex3a,
        )
data = np.load('convergence_heat.npz')
for name in data.files:
    exec(name + '=data[\''+name+'\']')
plt.figure(figsize=[6.4*1.7, 4*1.7])
plt.loglog(dts,err_be,'C1',label = 'BE')
plt.scatter(dts[False==change_be],err_be[False==change_be],color='C1',marker=plot_fkt.marker(-1))

plt.loglog(dts,err_ex3,'C2',label = 'BE 3 extrapolation')
plt.scatter(dts[False==change_ex3],err_ex3[False==change_ex3],color='C2',marker=plot_fkt.marker(-1))

plt.loglog(dts,err_ex3c,'C3',label = 'convex adaptation')
plot_fkt.plot_markers(dts,err_ex3c,False==change_ex3c,order_ex3c,'C3',label=False)

plt.loglog(dts,err_ex3a,'C4',label = 'free adaptation',zorder=2.1)
plot_fkt.plot_markers(dts,err_ex3a,False==change_ex3a,order_ex3a,'C4',label=False,zorder=2.1)

draw_order(x0=1.2e-5,y0=1e-3,x1=1.2e-4,p=1,ftext=1.2,text='above')
draw_order(x0=2e-5,y0=0.6e-9,x1=2e-4,p=3,ftext=1,text='below')
draw_order(x0=0.6e-4,y0=2e-6,x1=0.6e-3,p=2.3,ftext=1.2,text='above')

show_legend_marker(order=[3,2,1], loc='lower right')
plt.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
plt.grid()

plt.xlabel('$\Delta t$')
plt.ylabel('$e_{rel}$')

plt.savefig('conv_heat.pdf',bbox_inches = "tight")
np.min(dts[change_ex3a==True])<jupyter_output><empty_output><jupyter_text>Execution time<jupyter_code>np.max(dts[np.logical_and(change_ex3a==False, dts<1e-4)])
solver = Solver(rkm = ex3,
               dt = dt,
               t_final = 0.01,
               b_fixed=False,
               tol_neg=1e-8,
               tol_change = np.inf,
               p = [3,2,1],
               theta = [1],
               solver = cp.MOSEK,
               solver_eqs=solver_Matrix,
               convex = False,
               LP_opts = {'reduce':True})

problem_ADR = Problem(f=A_heat,
                 u0=u0,
                 minval=0,
                 maxval=np.inf)

t_adap = timeit.timeit(lambda:RK_integrate(solver=solver,problem=problem_ADR,verbose=False),number=30)

solver = Solver(rkm = ex3,
               dt = 2.9763514416313192e-05,
               t_final = 0.01,
               b_fixed=False,
               tol_neg=1e-8,
               tol_change = np.inf,
               p = [3,2,1],
               theta = [1],
               solver = cp.MOSEK,
               solver_eqs=solver_Matrix,
               convex = False,
               LP_opts = {'reduce':True})

t_small_dt = timeit.timeit(lambda:RK_integrate(solver=solver,problem=problem_ADR,verbose=False),number=30)

display('adaption')
display(t_adap)

display('small dt')
display(t_small_dt)<jupyter_output><empty_output><jupyter_text>Advection-decay equation<jupyter_code>def f_upwind_advection_influx(t,u):
    dx = 1/len(u)
    du = np.zeros_like(u)
    du[1:] = -(u[1:]-u[:-1])
    du[0] = -u[0]+1 #boundary condition
    return du/(dx)

def f_adv_dec(t,u):
    a = 1 #advection constant
    K = 1 #decay constant K =10
    T_a = 0
    
    du = a *f_upwind_advection_influx(t,u)+K*(T_a-u)
    return du

N = 100
dx = 1/N
dt = 1.7*dx
dt = 1.5*dx
print(dt)

solver = Solver(rkm = dp5,
               dt = dt,
               t_final = 1,
               b_fixed=False,
               tol_neg=1e-10,
               tol_change = 5,
               p = [4,3,2],
               theta = [1],
               solver = cp.MOSEK,
               LP_opts = {'reduce':True})

u0 = np.zeros(N)
#u0[0]=1
problem_adv_dec = Problem(f=f_adv_dec,
                 u0=u0,
                 minval=0,
                 maxval=np.inf)

status,t,u,b = RK_integrate(solver=solver,problem=problem_adv_dec,verbose=True)

t = np.array(t)
u = np.array(u).T
b = np.array(b).T
utils.show_status(status)
index = [0,13,27,40,53,-1]
x=np.linspace(0,1,100)
plt.figure(figsize=[6.4, 4])
for i in index:
    plt.plot(x,u[:,i],label='t='+str(round(t[i],2)))
    
plt.grid()
plt.xlabel('$x$')
plt.ylabel('$u$')
plt.legend(loc="center left", bbox_to_anchor=(1.0, 0.5))
#plt.title('Advection-Decay $\Delta t ='+str(dt)+'$')
plt.savefig('Advection_Decay.pdf',bbox_inches = "tight")
b_org = dp5.b.copy()
b_org.shape=(len(b_org),1)

order = np.array(status['order'])

plt.figure(figsize=[6.4, 2])
norm = np.linalg.norm((b-b_org),axis=0,ord=1)
plt.plot(t[1:],norm[1:],'-',color='C1')
mask=(norm[1:]==0)
plot_fkt.plot_markers(t[1:],norm[1:],mask,order[1:],'C1')
plt.xlabel('$t$')
plt.ylabel(r'$\| \tilde{b} - b \|_1 $')
plt.grid()
plt.legend(loc="center left", bbox_to_anchor=(1.0, 0.5))
plt.savefig('b_Advection_Decay.pdf',bbox_inches = "tight")<jupyter_output><empty_output><jupyter_text>Convergence<jupyter_code>def ref_adv_dec(N,t,dt=None,a=1,K=1):
    """
    Calculatest the exact solutio using a matrix exponential
    if no dt is given returns only the final value
    
    for more details on how the solution is constructed see seperate notebook
    """
    
    dx = 1/N

    A_adv = 1/dx * (-np.diag(np.ones(N))+np.diag(np.ones(N-1),-1))
    A_adv_dec = a*A_adv - K*np.diag(np.ones(N))
    L = np.zeros([N+1,N+1])
    L[1:,1:] = A_adv_dec #New right hand side
    L[1,0] = a/dx

    mu = np.zeros(N+1)#u_0 and initial Point
    mu[0] = 1
    

    if dt:
        t = np.arange(0,t+dt,dt)
        u = np.zeros([N,len(t)])
        u[:,0] = np.zeros(N)
        for i in range(1,len(t)):
            t_hat = t[i]
            u_hat = linalg.expm(t_hat*L)@mu
            u[:,i] = u_hat[1:]
        return t,u

    
    else:
        u_hat = linalg.expm(t*L)@mu
        return t,u_hat[1:]
    
  

N = 100
dx = 1/N
dt = 1.7*dx

dts = np.logspace(-2.2,np.log10(dt),num=40)
#dts = dts[dts<dt]
t_end = 0.5


#reference

reference = np.zeros([N,len(dts)])
for i in range(len(dts)):
    reference[:,i]= ref_adv_dec(N,t_end)[1]


u0 = np.zeros(N)
#u0[0]=1
problem_adv_dec = Problem(f=f_adv_dec,
                 u0=u0,
                 minval=0,
                 maxval=np.inf)

solver_fix = Solver(rkm = dp5,
               dt = dt,
               t_final = t_end,
               b_fixed=True,
               tol_neg=1e-11,
               tol_change = 5,
               p = [4,3,2,1],
               theta = [1],
               solver = cp.MOSEK,
               LP_opts = {})

solver_adp = Solver(rkm = dp5,
               dt = dt,
               t_final = t_end,
               b_fixed=False,
               tol_neg=1e-11,
               tol_change = np.inf,
               p = [4,3,2],
               theta = [1],
               solver = cp.MOSEK,
               LP_opts = {})

print('fixed')
sol_fix,err_fix,change_fix = utils.plot_convergence(problem_adv_dec,solver_fix,dts,reference,error='rel',step = -1,
                              Params=dict(),get_order=False)

print('adapt')
sol_adp,err_adp,change_adp,order_adp = utils.plot_convergence(problem_adv_dec,solver_adp,dts,reference,error='rel',step = -1,
                              Params=dict(),get_order=True)
np.savez('convergence_adde.npz', 
         dts=dts,
        sol_fix = sol_fix,err_fix = err_fix,change_fix = change_fix,
        sol_adp = sol_adp,err_adp=err_adp,change_adp=change_adp,order_adp=order_adp)
data = np.load('convergence_adde.npz')
for name in data.files:
    exec(name + '=data[\''+name+'\']')
plt.figure(figsize=[6.4*1.7, 4*1.7])
plt.loglog(dts,err_fix,'C0',label = 'Dormand--Prince RK5',zorder=1.9)
plt.loglog(dts,err_adp,'C1',label = 'free adaptation',zorder=2)



plt.scatter(dts[False==change_fix],err_fix[False==change_fix],zorder=2,marker=plot_fkt.marker(-1))
plot_fkt.plot_markers(dts,err_adp,False==change_adp,order_adp,'C1',label=False)


draw_order(x0=0.8e-2,x1=1.5e-2,y0=0.15e-7,p=5,ftext=1.2,text='below')


show_legend_marker(order=[4,3], loc='center left')
plt.legend()


plt.grid(which='both')
ax = plt.gca()
ax.set_xticks((7e-3,8e-3,9e-3,1e-2,1.5e-2))
ax.set_xticklabels(('$7x10^{-3}$','$8x10^{-3}$','$9x10^{-3}$','$10^{-2}$','$1.5x10^{-2}$'))

plt.xlabel('$\Delta t$')
plt.ylabel('$e_{rel}$')


plt.savefig('conv_adde.pdf', bbox_inches="tight")
#Store $b$ to plot stability region
np.save('b_ex_adde.npy', b[:,np.nonzero(dp5.b[0]!=b[0,:])])
plt.plot(sol_adp[-2])
#plt.plot(reference[:,2])
print(np.min(dts[change_adp == True]))
print(np.max(dts[change_adp == False]))
display(dts[-2])<jupyter_output><empty_output><jupyter_text>Execution time<jupyter_code>dt = 1.5*dx
solver = Solver(rkm = dp5,
               dt = dt,
               t_final = 1,
               b_fixed=False,
               tol_neg=1e-10,
               tol_change = 5,
               p = [4,3,2],
               theta = [1],
               solver = cp.MOSEK,
               LP_opts = {'reduce':True})

u0 = np.zeros(N)
problem_adv_dec = Problem(f=f_adv_dec,
                 u0=u0,
                 minval=0,
                 maxval=np.inf)

t_adap = timeit.timeit(lambda:RK_integrate(solver=solver,problem=problem_adv_dec,verbose=False),number=30)

solver = Solver(rkm = dp5,
               dt = 8e-3,
               t_final = 1,
               b_fixed=False,
               tol_neg=1e-10,
               tol_change = 5,
               p = [4,3,2],
               theta = [1],
               solver = cp.MOSEK,
               LP_opts = {'reduce':True})

t_small_dt = timeit.timeit(lambda:RK_integrate(solver=solver,problem=problem_adv_dec,verbose=False),number=30)

display('adaption')
display(t_adap)

display('small dt')
display(t_small_dt)

#status,t,u,b = RK_integrate(solver=solver,problem=problem_adv_dec,verbose=True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17269599666388658, 'fraction_numerical': 0.038313177648040035, 'mean_word_length': 2.8781079442086113, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '277151', 'license': 'MIT', 'n_tokens_mistral': 8413, 'n_tokens_neox': 8012, 'n_words': 1120, 'path': 'Convergence.ipynb', 'repo_name': 'ketch/Positive_RK_Reproducibility', 'size': '19184'}