Step Size

Hans Scharler
/
Step Size
about 11 hours ago
2
13
6
0
9883
game.m
function g=game()
    % STEP SIZE: you are the step-size controller of an ODE solver.
    g.init=@init; g.update=@update; g.render=@render;
end
function s=init(ax)
    ax.XLim=[0 800]; ax.YLim=[0 600]; ax.YDir='normal'; ax.Color=[1 1 1];
    % Course: S-bends [position height width] and stiff zones [start end].
    s.bend=[1400 120 110;2000 -200 90;2600 160 70;3900 -150 55;4300 200 60;...
        4700 -180 55;5350 100 80;6000 -120 50;6350 150 50];
    s.stiff=[3000 3400;5100 5500]; s.L=7040; s.best=inf;
    gray=[.55 .58 .62]; blue=[0 .447 .741]; orange=[.85 .325 .098];
    wall=[.93 .94 .95]; dark=[.15 .15 .15]; mid=[.3 .3 .3];
    G.zone=poly(ax,[1 .93 .82],'none');
    G.field=ln(ax,[.66 .7 .76],.8,'-');
    G.exact=ln(ax,[.494 .184 .556],1.2,':');
    G.top=poly(ax,wall,gray);
    G.bot=poly(ax,wall,gray);
    G.label=[txt(ax,0,0,'STIFF',11,orange) txt(ax,0,0,'STIFF',11,orange)];
    G.finish=ln(ax,[.47 .67 .19],3,'-');
    G.reject=ln(ax,[.8 .1 .1],1.5,'none','x',9);
    G.trial=ln(ax,orange,2,'--','.',14);
    G.chord=ln(ax,blue,1.5,'-');
    G.trail=ln(ax,blue,1.5,'-','o',4);
    G.probe=ln(ax,blue,1,'none','o',11);
    G.note=txt(ax,400,512,'',12,[.8 .1 .1]);
    poly(ax,[1 1 1],'none',[0 800 800 0],[532 532 600 600]);
    txt(ax,22,580,'STEP SIZE',20,dark,'left');
    txt(ax,23,555,'you are the step-size controller',10,gray,'left');
    G.h=txt(ax,400,575,'',22,blue);
    G.info=txt(ax,400,548,'',10,mid);
    ln(ax,gray,.5,'-','none',1,[590 780 780 590 590],[544 544 592 592 544]);
    G.hist=ln(ax,blue,1.2,'-','none',1,590+4.75*floor((1:80)/2),zeros(1,80));
    txt(ax,585,568,'log2 h',8,gray,'right');
    poly(ax,[1 1 1],'none',[0 800 800 0],[0 0 30 30]);
    txt(ax,400,16,'SPACE / CLICK: reject the orange trial step (h halves)',10,mid);
    G.panel=poly(ax,[.98 .98 .98],gray,[90 710 710 90],[95 95 505 505]);
    G.title=txt(ax,400,455,'',24,dark);
    G.body=txt(ax,400,320,'',11,[.2 .2 .2]); G.body.FontName='Consolas';
    G.prompt=txt(ax,400,155,'',15,blue);
    s.g=G; s=reset(s); s.mode=0;
end
function h=poly(ax,c,edge,x,y)
    if nargin<4, x=nan; y=nan; end
    h=patch(ax,x,y,c,'EdgeColor',edge,'HitTest','off','PickableParts','none');
end
function h=ln(ax,c,w,style,mk,ms,x,y)
    if nargin<7, x=nan; y=nan; end
    if nargin<5, mk='none'; ms=6; end
    h=line(ax,x,y,'Color',c,'LineWidth',w,'LineStyle',style,'Marker',mk,...
        'MarkerSize',ms,'HitTest','off','PickableParts','none');
end
function h=txt(ax,x,y,str,sz,c,align)
    if nargin<7, align='center'; end
    h=text(ax,x,y,str,'Color',c,'FontSize',sz,'HorizontalAlignment',align,...
        'HitTest','off','PickableParts','none');
end
function y=center(s,x)
    % Exact solution: smooth S-bends plus a wiggle inside each stiff zone.
    b=s.bend; y=300+sum(b(:,2).*(1+tanh((x-b(:,1))./b(:,3)))/2,1);
    for k=1:2
        u=(x-s.stiff(k,1))/diff(s.stiff(k,:));
        y=y+(u>0&u<1).*8.*sin(pi*u).^2.*sin(2*pi*(x-s.stiff(k,1))/200);
    end
end
function w=halfWidth(x)
    w=interp1([-1e4 3000 5800 1e5],[46 40 34 34],x);
end
function m=slope(s,x,y)
    % The ODE: y' = c'(x) - k e, with e = y - c. The center line c solves it.
    % k is high in stiff zones and grows near the walls.
    e=y-center(s,x);
    k=.004+.03*any(x>s.stiff(:,1)&x<s.stiff(:,2),1)+.01*(e./halfWidth(x)).^2;
    m=center(s,x+.5)-center(s,x-.5)-k.*e;
end
function u=hitAt(s,x0,y0,m,len)
    % Sample after the start point, which was accepted by the previous step.
    n=ceil(len/3); u=len*(1:n)/n;
    i=find(abs(y0+m*u-center(s,x0+u))>halfWidth(x0+u)-6,1);
    if isempty(i), u=inf; else, u=u(i); end
end
function s=reset(s)
    s.mode=1; s.age=0; s.clock=0; s.h=8; s.tapped=false; s.flash=0; s.freeze=0;
    s.steps=0; s.failed=0; s.evals=0; s.strikes=0; s.reason={};
    s.tx=nan(1,30); s.ty=nan(1,30); s.rx=nan(1,8); s.ry=nan(1,8);
    s.hist=-4*ones(1,40); s.tx(end)=0; s.ty(end)=center(s,0);
    s=commit(s,0,center(s,0));
end
function s=commit(s,x0,y0)
    % Take the trial step; the next trial doubles unless this one was rejected.
    s.x0=x0; s.y0=y0; s.m=slope(s,x0,y0); s.len=s.h; s.u=0;
    s.speed=s.h/(.55+s.h/900); s.hit=hitAt(s,x0,y0,s.m,s.len);
    s.xe=x0+s.len; s.ye=y0+s.m*s.len; s.mt=slope(s,s.xe,s.ye);
    if ~s.tapped, s.h=min(256,2*s.h); end
    s.tapped=false; s.evals=s.evals+1;
end
function s=update(s,dt,pressed)
    s.age=s.age+dt;
    if s.mode~=1
        if pressed && (s.mode==0 || s.age>.7), s=reset(s); end
        return
    end
    s.clock=s.clock+dt; s.flash=max(0,s.flash-dt);
    if s.clock>=90
        s=stop(s,2,{sprintf('Warning: Wall-clock limit reached at t=%e.',(s.x0+s.u)/128)}); return
    end
    if pressed && s.h>8
        % Reject the trial step: same direction, half the length.
        s.failed=s.failed+1; s.evals=s.evals+1;
        s.rx=[s.rx(2:end) s.xe+s.h]; s.ry=[s.ry(2:end) s.ye+s.mt*s.h];
        s.h=s.h/2; s.tapped=true;
    end
    if s.freeze>0, s.freeze=s.freeze-dt; return, end
    s.u=s.u+s.speed*dt;
    if s.u>=s.hit
        s.strikes=s.strikes+1; s.failed=s.failed+1; s.flash=.8;
        s.rx=[s.rx(2:end) s.x0+s.hit]; s.ry=[s.ry(2:end) s.y0+s.m*s.hit];
        if s.strikes==3
            s=stop(s,2,{sprintf('Warning: Failure at t=%e.  Unable to meet',(s.x0+s.hit)/128);...
                'integration tolerances without reducing the step size';...
                'below the smallest value allowed (0.0625) at time t.'}); return
        end
        s.h=max(8,s.len/4); s.tapped=true; s.freeze=.6; s=commit(s,s.x0,s.y0); return
    end
    if s.u>=s.len
        s.steps=s.steps+1; s.tx=[s.tx(2:end) s.xe]; s.ty=[s.ty(2:end) s.ye];
        s.hist=[s.hist(2:end) log2(s.len/128)];
        if s.xe>=s.L, s.best=min(s.best,s.clock); s=stop(s,3,{}); return, end
        s=commit(s,s.xe,s.ye);
    end
end
function s=stop(s,mode,reason)
    s.mode=mode; s.age=0; s.reason=reason;
end
function render(~,s)
    G=s.g; u=min(s.u,s.len); px=s.x0+u; py=s.y0+s.m*u; cam=px-220;
    xs=cam-15+(0:80)*10.4; c=center(s,xs); w=halfWidth(xs);
    set(G.top,'XData',[xs-cam 840 -40],'YData',[c+w 600 600]);
    set(G.bot,'XData',[xs-cam 840 -40],'YData',[c-w 0 0]);
    set(G.exact,'XData',xs-cam,'YData',c);
    set(G.zone,'XData',s.stiff(:,[1 2 2 1])'-cam,'YData',[0;0;600;600]+[0 0]);
    for k=1:2
        a=mean(s.stiff(k,:)); G.label(k).Position=[a-cam center(s,a)+75 0];
    end
    % Slope field on a grid that scrolls with the camera.
    [gx,gy]=meshgrid(cam-mod(cam,40)+(0:21)*40,50:40:530);
    gx=gx(:)'; gy=gy(:)'; m=slope(s,gx,gy); n=9./sqrt(1+m.^2);
    set(G.field,'XData',reshape([gx-n;gx+n;nan*gx],1,[])-cam,...
        'YData',reshape([gy-n.*m;gy+n.*m;nan*gy],1,[]));
    set(G.finish,'XData',[1 1]*(s.L-cam),'YData',center(s,s.L)+[-50 50]);
    set(G.trail,'XData',[s.tx px]-cam,'YData',[s.ty py]);
    set(G.chord,'XData',[s.x0 s.xe]-cam,'YData',[s.y0 s.ye]);
    set(G.trial,'XData',[s.xe s.xe+s.h]-cam,'YData',[s.ye s.ye+s.mt*s.h]);
    set(G.reject,'XData',s.rx-cam,'YData',s.ry);
    c=[0 .447 .741]; if s.flash>0, c=[.8 .1 .1]; end
    set(G.probe,'XData',px-cam,'YData',py,'Color',c,'MarkerFaceColor',c);
    G.hist.YData=548+(repelem(s.hist,2)+4)*8;
    G.h.String=sprintf('h = %g',s.h/128);
    G.info.String=sprintf('clock %4.1f s  |  t = %4.1f / %d  |  failed %d  |  strikes %d / 3',...
        s.clock,min(px,s.L)/128,s.L/128,s.failed,s.strikes);
    G.note.String=''; if s.flash>0, G.note.String='STEP FAILED: hit the wall, h / 4'; end
    vis='off'; if s.mode~=1, vis='on'; end
    set([G.panel G.title G.body G.prompt],'Visible',vis);
    G.trial.Visible=s.mode==1;
    stats={sprintf('%d successful steps',s.steps);sprintf('%d failed attempts',s.failed);...
        sprintf('%d function evaluations',s.evals)};
    if s.mode==0
        G.title.String='STEP SIZE';
        G.body.String={'A probe follows the slope field in straight Euler steps.';'';...
            'Orange is the next trial step. Clean steps double h.';...
            'TAP to reject the trial: h halves. Tap again to halve again.';'';...
            'Big steps are fast. Too big, and they cut into the wall.';...
            'Shaded zones are STIFF: big steps overshoot and blow up.';'';...
            'Reach t = 55 within 90 s of wall-clock time.'};
        G.prompt.String='TAP TO INTEGRATE';
    elseif s.mode==2
        G.title.String='INTEGRATION FAILED';
        G.body.String=[s.reason;{''};stats];
        G.prompt.String='TAP TO TRY AGAIN';
    elseif s.mode==3
        G.title.String='INTEGRATION COMPLETE';
        G.body.String=[{sprintf('Reached t = %d in %.1f s of wall-clock time.',s.L/128,s.clock);...
            sprintf('Best this session: %.1f s',s.best);''};stats];
        G.prompt.String='TAP TO INTEGRATE AGAIN';
    end
end
Play Game
Preview

Play Game opens this in MATLAB Drive. Click Run this app there to play.