% Examination of the behavior of a second-order linear constant
% coefficient equation, looking at the characteristic equation,
% its roots, and the graph of the solution to the differential
% equation.  This demo animates the variation of one parameter to
% show the effect of this on the characteristic equation, the
% roots, and the solution.
%
% by Justin Dunmyre
% rewritten by Gavin LaRose, 1/2/2013

% (c)2012 Regents of the University of Michigan
% This work is licensed under the
% Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% or send a letter to
% Creative Commons,
% 444 Castro Street, Suite 900,
% Mountain View, CA, 94041, USA.

function Char_Eqn_Animated()
%-------------------------------------------------------------------------------
% configuration variables

% our equation is 
%    m x'' + c x' + k x = 0,
% with initial conditions
%    x(0) = x0, x'(0) = xp0
% good values may be m = 1, c = 0.5, k = 3; and m = 1, c = 3, k = 6
    m = 5;
    c = 22;
    k = 45;
    x0 = 2;
    xp0 = 0;

% t values on which to graph the solutions
    t_range = [0 10];
% default y-range on which to graph the solutions
    y_range = [ -5 5 ];
    
% we vary one of the parameters, m, c or k:
    param_to_vary = 'k';
% range of values to consider
    param_range = [ 5 45 ];
% number of steps to consider between these values
    num_param_values  = 50;

% color choices    
    real_color   = 'b';
    repeat_color = 'r';
    imag_color   = 'g';
    bggray = [.6 .6 .6];

%-------------------------------------------------------------------------------
% initialization
    
% initial values of the parameters
    startm = m;
    startc = c;
    startk = k;

% we reset one of these based on what we're varying
    if param_to_vary == 'm'
        startm = param_range(1)
        m = param_range(1)
        char_roots2 = roots( [ param_range(2) c k ] );
    elseif param_to_vary == 'c'
        startc = param_range(1)
        c = param_range(1)
        char_roots2 = roots( [ m param_range(2) k ] );
    else
        startk = param_range(1);
        k = param_range(1);
        char_roots2 = roots( [ m c param_range(2) ] );
    end
    char_roots1 = roots( [ m c k ] );
% the values we use for the parameter we're varying
    param_values = linspace( param_range(1), param_range(2), num_param_values );

% figure out reasonable ranges for axes
%     for the polynomial graph
    poly_x_min = min( vertcat(real(char_roots1), real(char_roots2)) ) - 1;
    if poly_x_min >= 0
        poly_x_min = -1;
    else
        poly_x_min = poly_x_min - 1;
    end
    poly_x_max = max( vertcat(real(char_roots1), real(char_roots2)) ) + 1;
    if poly_x_max <= 0
        poly_x_max = 1;
    else
        poly_x_max = poly_x_max + 1;
    end
    poly_x_range = [ poly_x_min poly_x_max ];

    x_mid = (poly_x_range(1) + poly_x_range(2))/2;
    poly_y_min = round(m*x_mid*x_mid + c*x_mid + k) - 2;
    tmpval1 = m*poly_x_range(1)*poly_x_range(1) + c*poly_x_range(1) + k;
    tmpval2 = m*poly_x_range(2)*poly_x_range(2) + c*poly_x_range(2) + k;
    evalc([param_to_vary '=' num2str(param_range(2))]);
    tmpval3 = m*poly_x_range(1)*poly_x_range(1) + c*poly_x_range(1) + k;
    tmpval4 = m*poly_x_range(2)*poly_x_range(2) + c*poly_x_range(2) + k;
    evalc([param_to_vary '=' num2str(param_range(1))]);
    poly_y_max   = max([tmpval1 tmpval2 tmpval3 tmpval4]) + 5;
    poly_y_range = [ poly_y_min poly_y_max ];
    
%    for the roots graph
    rvals1 = sort(char_roots1);
    rvals2 = sort(char_roots2);
    root_x_range = [ min( [ real(rvals1(1)), real(rvals2(1)) ] ) - 1 ...
                     max( [ real(rvals1(2)), real(rvals2(2)) ] ) + 1 ];
    root_y_range = [ min( [ imag(rvals1(1)), imag(rvals2(1)) ] ) - 1 ...
                     max( [ imag(rvals1(2)), imag(rvals2(2)) ] ) + 1 ];
    
%-------------------------------------------------------------------------------
% set up the graph

    scrsz = get(0,'ScreenSize');
    figh = figure('Position',[1 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2]);
    axes('Parent', figh, 'Color', bggray);
    hold on;

% set up the subplots
    poly_sub_plot = subplot(2,2,1, 'Parent', figh, 'Color', bggray);
    root_sub_plot = subplot(2,2,2, 'Parent', figh, 'Color', bggray);
    soln_sub_plot = subplot(2,2,[3 4], 'Parent', figh, 'Color', bggray);
    hold( poly_sub_plot );
    hold( root_sub_plot );
    hold( soln_sub_plot );

    axis( poly_sub_plot, [ poly_x_range poly_y_range ] );
    axis( root_sub_plot, [ root_x_range root_y_range ] );
    axis( soln_sub_plot, [ t_range y_range ] );
    
% add dashed lines to show axes
    plot( poly_x_range, [0 0], 'w:', 'Parent', poly_sub_plot );
    plot( [0 0], poly_y_range, 'w:', 'Parent', poly_sub_plot );
    plot( root_x_range, [0 0], 'w:', 'Parent', root_sub_plot );
    plot( [0 0], root_y_range, 'w:', 'Parent', root_sub_plot );
    plot( t_range, [0 0], 'w:', 'Parent', soln_sub_plot );
    
% add some titles and axis labels
    poly_title = title( ['Characteristic Eqn: 0 = ', num2str(m), ' r^2 + ' ...
                         num2str(c), ' r + k'], 'FontSize', 16, ...
                        'Parent', poly_sub_plot );
    poly_xlabel = xlabel( 'r', 'FontSize', 16, 'Parent', poly_sub_plot );
    poly_ylabel = ylabel( 'y', 'FontSize', 16, 'Parent', poly_sub_plot );
    root_title = title( 'Roots of Char Eqns', 'FontSize', 16, ...
                        'Parent', root_sub_plot );
    root_xlabel = xlabel( 'Re(r)', 'FontSize', 16, 'Parent', root_sub_plot );
    root_ylabel = ylabel( 'Im(r)', 'FontSize', 16, 'Parent', root_sub_plot );
    soln_title = title( ['Solution to ', num2str(m), 'x\prime\prime + ' ...
                         num2str(c), 'x\prime + k x = 0, x(0) = ', ...
                         num2str(x0), ', x\prime(0) = ', num2str(xp0)], ...
                        'Parent', soln_sub_plot, 'FontSize', 16 );
    soln_xlabel = xlabel( 't', 'FontSize', 16, 'Parent', soln_sub_plot );
    soln_ylabel = ylabel( 'x', 'FontSize', 16, 'Parent', soln_sub_plot );
     
%-------------------------------------------------------------------------------
% populate the figures

% values of r to plot
    r_vals = poly_x_range(1):.1:poly_x_range(2);
    
    if imag(char_roots1(1)) ~= 0
        use_color = imag_color;
    elseif char_roots1(1) == char_roots1(2)
        use_color = repeat_color;
    else
        use_color = real_color;
    end

    pp_h = plot( r_vals, m*r_vals.*r_vals + c*r_vals + k, ...
                 'LineStyle', '-', 'Color', use_color, 'LineWidth', 2, ...
                 'Parent', poly_sub_plot );
    legend( pp_h, [param_to_vary ' = ' ...
                   sprintf('%4.2f', param_range(1))]);

    rp_h = plot( real(char_roots1), imag(char_roots1), 'Marker', '.', ...
                 'Color', use_color, 'MarkerSize', 20, 'LineWidth', 2, ...
                 'LineStyle', 'none', 'Parent', root_sub_plot );
    legend( rp_h, [param_to_vary ' = ' ...
                   sprintf('%4.2f', param_range(1))]);
    
    [t_pl, x_pl] = ode45( @(t,x) [ x(2); (-k*x(1) - c*x(2))/m ], ...
                          t_range, [x0; xp0] );
    sp_h = plot( t_pl, x_pl(:,1), 'LineStyle', '-', 'Color', use_color, ...
                 'LineWidth', 2, 'Parent', soln_sub_plot );
    legend( sp_h, [param_to_vary ' = ' ...
                   sprintf('%4.2f', param_range(1))]);
    
    getframe();
    
%-------------------------------------------------------------------------------
% presolve for all parameter values

    y = zeros(length(param_values), length(r_vals));

    for j=1:length(param_values)
        
        % set the parameter value
        evalc([param_to_vary '=' num2str(param_values(j))]);
        
        % check what color the graphs should be
        proots = roots( [m c k] );
        if imag(proots(1)) ~= 0
            color_list{j} = imag_color;
        elseif proots(1) == proots(2)
            color_list{j} = repeat_color;
        else
            color_list{j} = real_color;
        end
        
        % the polynomial
        y(j,:) = m*r_vals.*r_vals + c*r_vals + k;
        
        [tempt,tempsol]=ode45( @(t,x) [ x(2); (-k*x(1) - c*x(2))/m ], ...
                               t_range, [x0; xp0] );
        t{j}=tempt;
        sol{j}=tempsol(:,1);
        
        if rem( j, 10 ) == 0
            disp(['Solving ODE: ' num2str(j) ' of ' ...
                  num2str(length(param_values))]);
        end
    end
    
    disp('Press space to start animation');
    pause

%-------------------------------------------------------------------------------
% animation

    for j=2:length(param_values)
        pVal = param_values(j);

        delete( pp_h );
        delete( rp_h );
        delete( sp_h );
        
        evalc([param_to_vary '=' num2str(pVal)]);

        pp_h = plot( r_vals, y(j,:), [color_list{j} '-'], 'LineWidth', 2, ...
                     'Parent', poly_sub_plot );
        legend( pp_h, [param_to_vary ' = ' ...
                       sprintf('%4.2f', pVal)]);

        proots = roots( [m c k] );
        rp_h = plot( real(proots), imag(proots), [color_list{j} '.'], ...
                     'MarkerSize', 20, 'LineWidth', 2, ...
                     'Parent', root_sub_plot );
        legend( rp_h, [param_to_vary ' = ' ...
                       sprintf('%4.2f', pVal)]);

        sp_h = plot( t{j},sol{j}, [color_list{j} '-'], 'LineWidth', 2, ...
                     'Parent', soln_sub_plot );
        legend( sp_h, [param_to_vary ' = ' ...
                       sprintf('%4.2f',pVal)]);

        getframe();
        
    end
    
end

% all done
%-------------------------------------------------------------------------------
