Hi Dr. Ha and everyone!

Anybody found a way to do manual integration with good results here?

I have this program (ick.cmp) that won't calculate the areas of a
rectangular or a triangular area when the deltas are medium small. Not
infinitessimal at all.

If you run the test with the global variable named "test" set to 1 it
will miscalculate like crazy. When you set "test" to 2 it uses enormous
deltas but comes up with much more accurate values. I tried messing
with tol#=1e-6 and it's still pretty bad.

I'm also tossing in a graph program (SoL.cmp for "speed of light") to
show the sharp knee in the mass/velocity curve near the speed of light
just for fun.

It really is practically a brick wall.

I'm using CMAP to generate graphics as I explore physics and math a bit.

/*

Integral of 1/x between x=0 and x=1 is infinity defined?

*/

/////////////////////////////////////////////////////////////
float test=1; // set test=1 or 2 for small deltas or large
/////////////////////////////////////////////////////////////

DrawTriangle(float x1, float y1, float x2, float y2, float x3, float y3)
{
plotline(x1, y1, x2, y2);
plotline(x2, y2, x3, y3);
plotline(x3, y3, x1, y1);
}


main()
{

float x, y, oldx=1.0, oldy=1.0;
float x0=1.0, y0=1.0; // option for large deltas of 1/2, 1/4...
float delta=.5; // for smaller and smaller steps
float thresh=delta; // threshold for display
float a1, a2; // area accumulators for upper trgl and lower rect

clearplot();
setrange(0, 0, 1, 10);
setop(b,15);
setop(c, 0);
plotlinearscale("y", 10);
plotlinearscale("x", 10);
for (x=1; x>.05; x=x-1/128)
{
y=1/x;
plotline(oldx, oldy, x, y);
if (test==1)
{
a1=a1+((x-oldx)*(y-oldy))/2; // triangular area above
a2=a2+(x-oldx)*oldy; // rectangular area below
}
if (x<=thresh)
{
if (test==2)
{
a1=(x-x0)*(y-y0)/2;
a2=(x-x0)*(y0);
}
setop(c,12);
setop(f,12);
setop(p,2);
plotrectangle(x, y0, x0, 0);
setop(c,13);
DrawTriangle(x0,y0,x,y,x,y0);
setop(c,0);

putxy(x, y, abs(a1));
putxy(x,1, abs(a2));

delta=delta/2;
thresh=thresh-delta;
a1=0;
a2=0;
x0=x;
y0=y;
}
oldx=x;
oldy=y;
}
putxy(.3, 5, "Smallish deltas result in gross errors when integrating
manually.");
putxy(.4, 4, "Compare run when test=2 (delta ~ .5, .25...) to test=1 (delta =
1/128).");
putxy(.5, 3, "Current test setting is");
putxy(.68, 3, test);
if (test==1) {putxy(.7,3, "ICK!");}
}

/*
Physics/Relativity

Display mass increase as a function of velocity.


f(v)=1/(1-v^2/c^2)
where:
v is velocity
c is speed of light

Uses drawline to improve resolution for curve with sharp knee.

*/

func(float v) // percent of SoL
{
float tmp, c=300e6;
tmp=1/(1-(v^2/100^2)^(1/2));
return(tmp);
}


DrawGrid()
{
float xeff;
for (x=0; x<=100; x=x+10)
{
plotline(x, 0, x, 1000);
plotline(0, x*10, 100, x*10);
xeff=x*10; // avoid 0*mass display, should be 1*mass
if(xeff==0)
{
xeff=1;
}
putxy(-8, 1100, "(times mass)");
putxy(40, -120, "(% Speed of Light)");
putxy(-8, 20+xeff, xeff);
putxy(101, 20+xeff, xeff);
putxy(x-1, -50, x);
}
}


DrawCurve()
{
float oldx=0, oldy=0, stepsize=1;
for (i=0;i<100;i=i+stepsize)
{
newy=func(i);
if (newy<=1000) // keep it on the graph
{
plotline(oldx, oldy, i, newy);
}
oldx=i;
oldy=newy;
// adjust stepsize for better resolution as curve turns up
if (i>90) {stepsize=0.1;}
if (i>99) {stepsize=0.01;}
}
}


LabelAPoint(float x, float y, float val)
{
putxy(x, y, "(");
putxy(x+1, y, val);
putxy(x+3+floor(log10(val)), y, ")");
plotarrow(x+2,y-50,x,0,'E');
}

LabelPoints() // show numerical values at 50% and 90% SoL
{
LabelAPoint(50,100,func(50));
LabelAPoint(90, 100, func(90));
}

main()
{
clearplot();
setrange(-20, -200, 120, 1200); // leave larger blank margins
setop(B,15); // background= white

setop(c, 4); // grid = red
DrawGrid();
LabelPoints();

setop(C,0); // curve = black
// plot(i, 0, 99.9, func(i)); // not good enough resolution
DrawCurve();
}