# [NTG-context] OT: looking for metapost/fun examples

David Arnold darnold at northcoast.com
Sun Mar 27 06:27:40 CEST 2005

```Gerben,

I am so happy I can contribute as a way of thanking you for all the
work that you've done for us. Thanks. Here, try this one, which uses an
RK4 routine. I just tried it out in Texshop, so I know they both
compile.

%This file creates two figures associated with the
%system x'=f(x,y), y'=g(x,y)
%1. Plots the graphs of x(t) and y(t)
%2. Plots the graph of (x(t),y(t)) in the phase plane.

%verbatimtex
%\input mtplain
%etex

%Generate standard eps
prologues:=2;

beginfig(0);

%Place RHS of x'=f(t,x,y) here
def fxy(expr t, x, y)=
(0.4-0.01*y)*x
enddef;

%Place RHS of y'=g(t,x,y) here
def gxy(expr t, x, y)=
(-0.3+0.005*x)*y
enddef;

%Declare some variables
path q, trajx, trajy;
pair L, R, B, T, xt, yt;
numeric sx[], sy[];

%Initialize clipping window
a:=0; b:=40;   %left and right of viewing rectangle
c:=0; d:=150;  %bottom and top of viewing rectangle

%Initialize timespan
tstart:=a;
tstop:=b;

%Initialize number of points to be plotted
N:=500;

%Calculate time increment dt for Euler's method
dt:=(tstop-tstart)/N;

%Scaling factors for horizontal and vertical axes. Note that this
produces
%an image that is 2 inches by 2 inches.
(b-a)*ux=1.75in;
(d-c)*uy=1.75in;

%Clipping boundary
q=(a,c)--(b,c)--(b,d)--(a,d)--cycle;

%Use Runge-Kutta4 to create path (t,x(t))

%Choose initial condition
t:=tstart;
x:=40;
y:=20;
trajx:=(t,x);
forever:
sx1:=fxy(t,x,y);
sy1:=gxy(t,x,y);
sx2:=fxy((t+dt/2),(x+dt*sx1/2),(y+dt*sy1/2));
sy2:=gxy((t+dt/2),(x+dt*sx1/2),(y+dt*sy1/2));
sx3:=fxy((t+dt/2),(x+dt*sx2/2),(y+dt*sy2/2));
sy3:=gxy((t+dt/2),(x+dt*sx2/2),(y+dt*sy2/2));
sx4:=fxy((t+dt),(x+dt*sx3),(y+dt*sy3));
sy4:=gxy((t+dt),(x+dt*sx3),(y+dt*sy3));
x:=x+dt*(sx1+2*sx2+2*sx3+sx4)/6;
y:=y+dt*(sy1+2*sy2+2*sy3+sy4)/6;
t:=t+dt;
trajx:=trajx..(t,x);
exitif ((t>tstop) or (t>b) or (x<c) or (x>d));
endfor;

%Use Runge-Kutta4 to create path (t,y(t))

%Choose initial condition
t:=tstart;
x:=40;
y:=20;
trajy:=(t,y);
forever:
sx1:=fxy(t,x,y);
sy1:=gxy(t,x,y);
sx2:=fxy((t+dt/2),(x+dt*sx1/2),(y+dt*sy1/2));
sy2:=gxy((t+dt/2),(x+dt*sx1/2),(y+dt*sy1/2));
sx3:=fxy((t+dt/2),(x+dt*sx2/2),(y+dt*sy2/2));
sy3:=gxy((t+dt/2),(x+dt*sx2/2),(y+dt*sy2/2));
sx4:=fxy((t+dt),(x+dt*sx3),(y+dt*sy3));
sy4:=gxy((t+dt),(x+dt*sx3),(y+dt*sy3));
x:=x+dt*(sx1+2*sx2+2*sx3+sx4)/6;
y:=y+dt*(sy1+2*sy2+2*sy3+sy4)/6;
t:=t+dt;
trajy:=trajy..(t,y);
exitif ((t>tstop) or (t>b) or (y<c) or (y>d));
endfor;

%Draw the paths x(t) and y(t) and clip them to bounding box
draw trajx xscaled ux yscaled uy withcolor red;
draw trajy xscaled ux yscaled uy withcolor red dashed evenly;
clip currentpicture to (q xscaled ux yscaled uy);

%Label graph x(t) and initial condition
len:= 0.65*(length trajx);
xt:=point len of trajx;
label.urt(btex \$\scriptstyle x(t)\$ etex, (xt xscaled ux yscaled uy));

%Label graph y(t) and initial condition
len:= 0.5*(length trajy);
yt:=point len of trajy;
label.lrt(btex \$\scriptstyle y(t)\$ etex, (yt xscaled ux yscaled uy));

%Initialize left and right endpoints on time-axis
L=(a*ux,0);R=(b*ux,0);

%Draw and label t-axis
drawarrow L--R;
label.rt(btex \$\scriptstyle t\$ etex,(b*ux,0));

%Initialize bottom and top endpoints on time-axis
B=(0,c*uy);T=(0,d*uy);

%Draw and label vertical axis
drawarrow B--T;
label.lft(btex \$\scriptstyle 0\$ etex, B);
label.lft(btex \$\scriptstyle 150\$ etex, T);

endfig;

beginfig(2);

%Make some variables local
save ux, uy;

%Place RHS of x'=f(t,x,y) here
def fxy(expr t, x, y)=
(0.4-0.01*y)*x
enddef;

%Place RHS of y'=g(t,x,y) here
def gxy(expr t, x, y)=
(-0.3+0.005*x)*y
enddef;

%Declare some variables
path q, trajxy;
pair L, R, B, T;

%Initialize clipping window
a:=0; b:=150;   %left and right of viewing rectangle
c:=0; d:=100;  %bottom and top of viewing rectangle

%Initialize timespan
tstart:=a;
tstop:=b;

%Initialize number of points to be plotted
N:=500;

%Calculate time increment dt for Euler's method
dt:=(tstop-tstart)/N;

%Scaling factors for horizontal and vertical axes. Note that this
produces
%an image that is 2 inches by 2 inches.
(b-a)*ux=1.75in;
(d-c)*uy=1.75in;

%Clipping boundary
q=(a,c)--(b,c)--(b,d)--(a,d)--cycle;

%Use Runge-Kutta4 to create path (x(t),y(t))

%Choose initial condition
t:=tstart;
x:=40;
y:=20;
trajxy:=(x,y);
forever:
sx1:=fxy(t,x,y);
sy1:=gxy(t,x,y);
sx2:=fxy((t+dt/2),(x+dt*sx1/2),(y+dt*sy1/2));
sy2:=gxy((t+dt/2),(x+dt*sx1/2),(y+dt*sy1/2));
sx3:=fxy((t+dt/2),(x+dt*sx2/2),(y+dt*sy2/2));
sy3:=gxy((t+dt/2),(x+dt*sx2/2),(y+dt*sy2/2));
sx4:=fxy((t+dt),(x+dt*sx3),(y+dt*sy3));
sy4:=gxy((t+dt),(x+dt*sx3),(y+dt*sy3));
x:=x+dt*(sx1+2*sx2+2*sx3+sx4)/6;
y:=y+dt*(sy1+2*sy2+2*sy3+sy4)/6;
t:=t+dt;
trajxy:=trajxy..(x,y);
exitif ((t>tstop) or (t>b) or (x<a) or (x>b) or (y<c) or (y>d));
endfor;

%Draw the paths x(t) and y(t) and clip them to bounding box
draw trajxy xscaled ux yscaled uy withcolor red;
clip currentpicture to (q xscaled ux yscaled uy);

%Initialize left and right endpoints on x-axis
L=(a*ux,0);R=(b*ux,0);

%Draw and label x-axis
drawarrow L--R;
label.rt(btex \$\scriptstyle x\$ etex,(b*ux,0));
label.bot(btex \$\scriptstyle 0\$ etex,L);
label.bot(btex \$\scriptstyle 150\$ etex,R);

%Initialize bottom and top endpoints on y-axis
B=(0,c*uy);T=(0,d*uy);

%Draw and label vertical axis
drawarrow B--T;
label.rt(btex \$\scriptstyle y\$ etex,(0,d*uy));
label.lft(btex \$\scriptstyle 0\$ etex, B);
label.lft(btex \$\scriptstyle 100\$ etex, T);

endfig;

end;

On Mar 26, 2005, at 3:19 PM, Gerben Wierda wrote:

> I am trying to learn metapost/fun, inline in ConTeXt source. Some
> basic things are clear, but now the issue is metapost itself.
>
> For instance, I would like to plot a Fourier approximation of a block
> function.
>
> For instance, I would like to plot a gaussian spread.
>
> I am looking for examples on how to do this. I need to do a bit of
> programming here and these are my initial projects.
>