human
I am working on the file "SC1v2/11_ODE.ipynb".
The first task is:
So, how do we use the stability region for each method?
The key point here is that we can't define $h$ a priori since we may not know the corresponding $\lambda$, thus, in principle, we can estimate $\lambda$ or just try the solver for different $h$.
Can you write Python program for it?
assistant
#Here is the Python program you requested.

def plot_log_eq_stability_region(y0=1.2, N=20, T=10, ODEsolver=eulerMethod):
    
    # How did w get this? Why do I need it?
    def f_1D_prime(y):
        return np.array(1-2*y)
    
    t0=0        
    if ODEsolver==eulerMethod:
        c1,c2,c3='b','k','k'
    elif ODEsolver==RK2:
        c1,c2,c3='k','b','k'
    else:
        c1,c2,c3='k','k','b'
        
    t_times, y_output = ODEsolver(t0,T,N,y0,f)

    fig = plt.figure(figsize=(L,L/3))
    ax1 = fig.add_subplot(131)
    ax2 = fig.add_subplot(132)
    ax3 = fig.add_subplot(133)

    ax1.axis('equal')
    ax1.grid(True)
    ax1.set_title('Numerical Approximation')
    ax1.axis([0, T, -0.5, 2])
    ax1.plot(t_times,y_output,'k-')
    ax1.plot(t_times,y_output,'b.')

    # The next code is to plot in logscale unstable numerical simulations
    h = (T-t0)/N
    y_all=h*f_1D_prime(y_output)
    y_pos=y_all>0
    y_neg=np.logical_not(y_pos)
    ax3.semilogy(t_times[y_pos],y_all[y_pos],'.r',ms=20)
    ax3.semilogy(t_times[y_neg],-y_all[y_neg],'.b',ms=20)
    ax3.grid(True)
    ax3.set_title(r'$\log10(|y_i|)$')

    zplot2(W1,ax2,line_color=c1)
    zplot2(W2,ax2,line_color=c2)
    zplot2(W4,ax2,line_color=c3)
    k_lambdah=h*f_1D_prime(y_output)
    ax2.plot(np.real(k_lambdah),np.imag(k_lambdah),'.r',ms=20)
    ax2.axis('equal')
    ax2.axis([-5, 2, -3.5, 3.5])
    ax2.grid(True)
    ax2.set_title('Stability Region and $k=\lambda\,h$')

interact(plot_log_eq_stability_region,y0=(-2,5,0.1),N=(2,100,1), T=(1,20,1), ODEsolver=radio_button_ODEsolvers)