Commit 731e5e49 authored by David Saxton's avatar David Saxton

Use 2nd degree Newton-Cotes formula (Simpson's rule) for numerical evaluation

of integrals.

svn path=/trunk/KDE/kdeedu/kmplot/; revision=528407
parent 1a9156ab
......@@ -16,7 +16,6 @@ TODO
* Draw tangents in trace mode (pipesmoker)
* Printing (pipesmoker)
- different paper sizes
* Use a better(faster!) algorithm for drawing integrals numeric. It's not urgent anymore since the implemention of Euler's method is better now.
* An export dialog where you can set the size and enable/disable monocrome.
* More printing options.
......@@ -27,6 +26,7 @@ IN PROGRESS
DONE
=========================
* Use a better(faster!) algorithm for drawing integrals numeric. It's not urgent anymore since the implemention of Euler's method is better now.
* Fix the unpolished lines?
* Popupmenu and tracemode for parametric- and polar functions
* Function edit Dialog
......@@ -81,4 +81,4 @@ DONE
* DCOP
* Make it possible to set a min OR max range value for functions.
* Checkbox in the function list to show/hide functions.
* Ability to set min and max values for the sliders.
\ No newline at end of file
* Ability to set min and max values for the sliders.
......@@ -244,19 +244,21 @@ void View::draw( QPaintDevice * dev, PlotMedium medium )
double View::value( Equation * eq, Function::PMode mode, double x )
{
double dx = (m_xmax-m_xmin)/1e3;
switch ( mode )
{
case Function::Derivative0:
return XParser::self()->fkt( eq, x );
case Function::Derivative1:
return XParser::self()->a1fkt( eq, x, (m_xmax-m_xmin)/1e3 );
return XParser::self()->derivative1( eq, x, dx );
case Function::Derivative2:
return XParser::self()->a2fkt( eq, x, (m_xmax-m_xmin)/1e3 );
return XParser::self()->derivative2( eq, x, dx );
case Function::Integral:
return XParser::self()->euler_method( x, eq );
return XParser::self()->integral( eq, x, dx );
}
kWarning() << k_funcinfo << "Unknown mode!\n";
......@@ -407,6 +409,8 @@ void View::plotfkt(Function *ufkt, QPainter *pDC)
}
mflg=2;
x = dmin;
dx = base_dx;
if ( p_mode == Function::Integral )
{
if ( ufkt->integral_use_precision )
......@@ -415,28 +419,10 @@ void View::plotfkt(Function *ufkt, QPainter *pDC)
// else
// dx = ufkt->integral_precision;
startProgressBar((int)double((dmax-dmin)/dx)/2);
x = ufkt->eq[0]->oldx = ufkt->startx.value(); //the initial x-point
ufkt->eq[0]->oldy = ufkt->starty.value();
ufkt->eq[0]->oldyprim = ufkt->integral_precision;
}
else
x=dmin;
bool forward_direction;
if (dmin<0 && dmax<0)
forward_direction = false;
else
forward_direction = true;
while ((x>=dmin && x<=dmax) || (p_mode == Function::Integral && x>=dmin && !forward_direction) || (p_mode == Function::Integral && x<=dmax && forward_direction))
while ( x>=dmin && x<=dmax )
{
if ( p_mode == Function::Integral && stop_calculating)
{
p_mode = Function::Derivative1;
x=dmax+1;
continue;
}
p2 = dgr.toPixel( realValue( ufkt, p_mode, x ) );
if ( p_mode == Function::Integral && (int(x*100)%2==0) )
......@@ -489,34 +475,14 @@ void View::plotfkt(Function *ufkt, QPainter *pDC)
}
mflg=0;
}
if ( p_mode == Function::Integral )
{
if ( forward_direction)
{
x=x+dx;
if (x>dmax && p_mode == Function::Integral )
{
forward_direction = false;
x = ufkt->eq[0]->oldx = ufkt->startx.value();
ufkt->eq[0]->oldy = ufkt->starty.value();
ufkt->eq[0]->oldyprim = ufkt->integral_precision;
mflg=2;
}
}
else
x=x-dx; // go backwards
}
if ( dxTooBig )
dx *= 0.5;
else
{
if ( dxTooBig )
dx *= 0.5;
else
{
if ( dxTooSmall )
dx *= 2.0;
x=x+dx;
}
if ( dxTooSmall )
dx *= 2.0;
x=x+dx;
}
}
}
......@@ -544,7 +510,7 @@ QPen View::penForPlot( Function *ufkt, Function::PMode p_mode, bool antialias )
pen.setCapStyle(Qt::RoundCap);
// pen.setStyle( Qt::DashLine );
double lineWidth_mm;
double lineWidth_mm = 0.0;
switch ( p_mode )
{
......@@ -728,7 +694,7 @@ bool View::root(double *x0, Equation *it)
do
{
double df = XParser::self()->a1fkt( it, *x0, (m_xmax-m_xmin)/1e3 );
double df = XParser::self()->derivative1( it, *x0, (m_xmax-m_xmin)/1e3 );
// kDebug() << "df1="<<df<<endl;
// if ( qAbs(df) < 1e-6 )
// df = 1e-6 * ((df < 0) ? -1 : 1);
......@@ -1019,7 +985,7 @@ void View::getPlotUnderMouse()
m_currentFunctionPlot = Function::Derivative0;
m_trace_x = 0.0;
int best_id;
int best_id = -1;
double best_distance = 1e30; // a nice large number
QPointF best_cspos;
......@@ -1733,22 +1699,14 @@ QPointF View::findMinMaxValue(Function *ufkt, Function::PMode p_mode, ExtremaTyp
else
dx = 1.0;
}
dmin = ufkt->eq[0]->oldx = ufkt->startx.value(); //the initial x-point
ufkt->eq[0]->oldy = ufkt->starty.value();
ufkt->eq[0]->oldyprim = ufkt->integral_precision;
}
else
dx = (dmax-dmin)/area.width();
x=dmin;
bool forward_direction;
if (dmin<0 && dmax<0)
forward_direction = false;
else
forward_direction = true;
while ((x>=dmin && x<=dmax) || (p_mode == Function::Integral && x>=dmin && !forward_direction) || (p_mode == Function::Integral && x<=dmax && forward_direction))
while ( (x>=dmin && x<=dmax) )
{
if ( p_mode == Function::Integral && stop_calculating)
{
......@@ -1764,7 +1722,7 @@ QPointF View::findMinMaxValue(Function *ufkt, Function::PMode p_mode, ExtremaTyp
kDebug() << "y " << y << endl;
if (x>=dmin && x<=dmax)
{
if ( start)
if ( start )
{
result_x = x;
result_y = y;
......@@ -1794,31 +1752,15 @@ QPointF View::findMinMaxValue(Function *ufkt, Function::PMode p_mode, ExtremaTyp
}
}
}
if (p_mode==Function::Integral)
{
if ( forward_direction)
{
x=x+dx;
if (x>dmax && p_mode== Function::Integral )
{
forward_direction = false;
x = ufkt->eq[0]->oldx = ufkt->startx.value();
ufkt->eq[0]->oldy = ufkt->starty.value();
ufkt->eq[0]->oldyprim = ufkt->integral_precision;
}
}
else
x=x-dx; // go backwards
}
else
x=x+dx;
x += dx;
}
return QPointF( result_x, result_y );
}
double View::getYValue( Function *ufkt, Function::PMode p_mode, double x, double y, const QString &str_parameter )
double View::getYValue( Function *ufkt, Function::PMode p_mode, double x, const QString &str_parameter )
{
assert( ufkt->type() == Function::Cartesian );
......@@ -1835,57 +1777,7 @@ double View::getYValue( Function *ufkt, Function::PMode p_mode, double x, double
}
}
if ( p_mode != Function::Integral )
return value( ufkt->eq[0], p_mode, x );
double dmin = getXmin( ufkt );
double dmax = getXmax( ufkt );
const double target = x; //this is the x-value the user had chosen
bool forward_direction;
if ( target>=0)
forward_direction = true;
else
forward_direction = false;
if(dmin==dmax) //no special plot range is specified. Use the screen border instead.
{
dmin=m_xmin;
dmax=m_xmax;
}
double dx;
if ( ufkt->integral_use_precision )
dx = ufkt->integral_precision*(dmax-dmin)/area.width();
else
dx=(dmax-dmin)/area.width();
x = ufkt->eq[0]->oldx = ufkt->startx.value(); //the initial x-point
ufkt->eq[0]->oldy = ufkt->starty.value();
ufkt->eq[0]->oldyprim = ufkt->integral_precision;
while ( x>=dmin )
{
y = XParser::self()->euler_method( x, ufkt->eq[0] );
if ( (x+dx > target && forward_direction) || ( x+dx < target && !forward_direction)) //right x-value is found
return y;
if (forward_direction)
{
x=x+dx;
if (x>dmax)
{
forward_direction = false;
x = ufkt->eq[0]->oldx = ufkt->startx.value();
ufkt->eq[0]->oldy = ufkt->starty.value();
ufkt->eq[0]->oldyprim = ufkt->integral_precision;
}
}
else
x=x-dx; // go backwards
}
kWarning() << k_funcinfo << "Could not find y value!\n";
return 0.0;
return value( ufkt->eq[0], p_mode, x );
}
void View::keyPressEvent( QKeyEvent * e )
......@@ -2064,17 +1956,8 @@ double View::areaUnderGraph( IntegralDrawSettings s )
}
double dx;
if ( s.pMode == Function::Integral )
{
if ( ufkt->integral_use_precision )
dx = ufkt->integral_precision*(s.dmax-s.dmin)/area.width();
else
dx = (s.dmax-s.dmin)/area.width();
s.dmin = ufkt->eq[0]->oldx = ufkt->startx.value(); //the initial x-point
ufkt->eq[0]->oldy = ufkt->starty.value();
ufkt->eq[0]->oldyprim = ufkt->integral_precision;
}
if ( (s.pMode == Function::Integral) && ufkt->integral_use_precision )
dx = ufkt->integral_precision*(s.dmax-s.dmin)/area.width();
else
dx = (s.dmax-s.dmin)/area.width();
......@@ -2084,12 +1967,9 @@ double View::areaUnderGraph( IntegralDrawSettings s )
int intervals = qRound( (s.dmax-s.dmin)/dx );
dx = (s.dmax-s.dmin) / intervals;
bool forward_direction = (s.dmin>=0) || (s.dmax>=0);
double calculated_area=0;
double x = s.dmin;
// while ( (x>=dmin && x<=dmax) || (p_mode == Function::Integral && x>=dmin && !forward_direction) || (p_mode == Function::Integral && x<=dmax && forward_direction))
// {
for ( int i = 0; i <= intervals; ++i )
{
double y = value( ufkt->eq[0], s.pMode, x );
......@@ -2100,24 +1980,7 @@ double View::areaUnderGraph( IntegralDrawSettings s )
else
calculated_area += dx*y;
if ( s.pMode == Function::Integral )
{
if ( forward_direction)
{
x=x+dx;
if (x>s.dmax && s.pMode == Function::Integral )
{
forward_direction = false;
x = ufkt->eq[0]->oldx = ufkt->startx.value();
ufkt->eq[0]->oldy = ufkt->starty.value();
ufkt->eq[0]->oldyprim = ufkt->integral_precision;
}
}
else
x=x-dx; // go backwards
}
else
x=x+dx;
x=x+dx;
}
m_integralDrawSettings = s;
......@@ -2279,7 +2142,7 @@ void View::invertColor(QColor &org, QColor &inv)
void View::updateCursor()
{
Cursor newCursor;
Cursor newCursor = m_prevCursor;
if ( isDrawing && (m_zoomMode != Translating) )
newCursor = CursorWait;
......
......@@ -114,7 +114,7 @@ class View : public QWidget, virtual public ViewIface
*/
QPointF findMinMaxValue( Function * function, Function::PMode p_mode, ExtremaType type, double dmin, double dmax,const QString & parameter );
/// get a y-value from a x-value
double getYValue( Function * function, Function::PMode p_mode, double x, double y, const QString & parameter );
double getYValue( Function * function, Function::PMode p_mode, double x, const QString & parameter );
/**
* Calculates the area between the given plot and the x-axis
* (from x = \p dmin to x = \p dmax). The area will also be colored in.
......
......@@ -29,6 +29,7 @@
#include <kdebug.h>
#include <assert.h>
#include <cmath>
......@@ -87,9 +88,6 @@ Equation::Equation( Type type, Function * parent )
{
mem = new unsigned char [MEMSIZE];
mptr = 0;
oldyprim = 0.0;
oldx = 0.0;
oldy = 0;
}
......@@ -241,8 +239,8 @@ Function::Function( Type type )
integral_precision = 1.0;
use_slider = -1;
startx.updateExpression( "0" );
starty.updateExpression( "0" );
m_startX.updateExpression( "0" );
m_startY.updateExpression( "0" );
// min/max stuff
dmin.updateExpression( QString("-")+QChar(960) );
......@@ -281,8 +279,8 @@ bool Function::copyFrom( const Function & function )
COPY_AND_CHECK( integral_use_precision ); // 4
COPY_AND_CHECK( dmin ); // 5
COPY_AND_CHECK( dmax ); // 6
COPY_AND_CHECK( starty ); // 7
COPY_AND_CHECK( startx ); // 8
COPY_AND_CHECK( m_startX ); // 7
COPY_AND_CHECK( m_startY ); // 8
COPY_AND_CHECK( integral_precision ); // 9
COPY_AND_CHECK( use_slider ); // 10
COPY_AND_CHECK( usecustomxmin ); // 11
......@@ -312,6 +310,16 @@ bool Function::copyFrom( const Function & function )
}
void Function::setIntegralStart( const Value & x, const Value & y )
{
assert( type() == Cartesian ); // Integral only applicable for cartesians
m_startX = x;
m_startY = y;
eq[0]->lastIntegralPoint = QPointF( 0.0, 0.0 );
}
QString Function::typeToString( Type type )
{
switch ( type )
......
......@@ -28,6 +28,7 @@
#define FUNCTION_H
#include <QColor>
#include <QPointF>
#include <QString>
......@@ -151,9 +152,9 @@ class Equation
*/
bool setFstr( const QString & fstr, bool force = false );
double oldx; ///< needed for Euler's method, the last x-value
double oldy; ///< The last y-value needed for Euler's method
double oldyprim; ///< needed for Euler's method, the last y'.value
QPointF lastIntegralPoint; ///< needed for numeric integration
double oldy;
protected:
const Type m_type;
......@@ -237,17 +238,17 @@ class Function
*/
Value dmax;
/**
* Start position for Euler's method, the initial y-value.
* Sets the start position and value of the integral.
*/
Value starty;
void setIntegralStart( const Value & x, const Value & y );
/**
* Start position for Euler's method, the initial x-value.
* The initial x-value in calculating integrals.
*/
Value startx;
Value integralInitialX() const { return m_startX; }
/**
* Precision when drawing numeric prime-functions.
* The initial y-value in calculating integrals.
*/
double integral_precision;
Value integralInitialY() const { return m_startY; }
/**
* -1: None (use list)
* else: slider number.
......@@ -259,10 +260,12 @@ class Function
QList<Value> parameters;
bool usecustomxmin:1;
bool usecustomxmax:1;
double integral_precision;
// TODO double slider_min, slider_max; ///< extreme values of the slider
protected:
const Type m_type;
Value m_startX, m_startY;
};
......
......@@ -347,8 +347,8 @@ void FunctionEditor::initFromCartesian()
m_editor->showIntegral->setChecked( f->integral.visible );
m_editor->customPrecision->setChecked( f->integral_use_precision );
m_editor->txtInitX->setText(f->startx.expression());
m_editor->txtInitY->setText(f->starty.expression());
m_editor->txtInitX->setText( f->integralInitialX().expression() );
m_editor->txtInitY->setText( f->integralInitialY().expression() );
m_editor->stackedWidget->setCurrentIndex( 0 );
m_editor->tabWidget->setCurrentIndex( 0 );
......@@ -554,13 +554,7 @@ void FunctionEditor::saveCartesian()
tempFunction.f0.color = m_editor->cartesian_f_lineColor->color();
tempFunction.integral.visible = m_editor->showIntegral->isChecked();
ok = tempFunction.startx.updateExpression( m_editor->txtInitX->text() );
if ( tempFunction.integral.visible && !ok )
return;
ok = tempFunction.starty.updateExpression( m_editor->txtInitY->text() );
if ( tempFunction.integral.visible && !ok )
return;
tempFunction.setIntegralStart( m_editor->txtInitX->text(), m_editor->txtInitY->text() );
tempFunction.integral.color = m_editor->cartesian_F_lineColor->color();
tempFunction.integral_use_precision = m_editor->customPrecision->isChecked();
......
......@@ -229,7 +229,8 @@ void KMinMax::cmdFind_clicked()
KMessageBox::sorry(this, i18n("Please choose a function"));
return;
}
double dmin, dmax;
double dmin = 0.0;
double dmax = 0.0;
dmin = XParser::self()->eval(m_mainWidget->min->text() );
if ( XParser::self()->parserError( true )!=0 )
{
......@@ -326,7 +327,7 @@ void KMinMax::cmdFind_clicked()
case CalculateY:
{
double value = View::self()->getYValue( ufkt, p_mode, dmin, dmax, parameter );
double value = View::self()->getYValue( ufkt, p_mode, dmin, parameter );
if ( !View::self()->isCalculationStopped() )
{
QString tmp;
......
......@@ -201,8 +201,8 @@ void KmPlotIO::addFunction( QDomDocument & doc, QDomElement & root, Function * f
tag.setAttribute( "integral-width", function->integral.lineWidth );
tag.setAttribute( "integral-use-precision", function->integral_use_precision );
tag.setAttribute( "integral-precision", function->integral_precision );
tag.setAttribute( "integral-startx", function->startx.expression() );
tag.setAttribute( "integral-starty", function->starty.expression() );
tag.setAttribute( "integral-startx", function->integralInitialX().expression() );
tag.setAttribute( "integral-starty", function->integralInitialY().expression() );
tag.setAttribute( "type", Function::typeToString( function->type() ) );
for ( unsigned i=0; i<2; ++i )
......@@ -517,8 +517,7 @@ void KmPlotIO::parseFunction( const QDomElement & n, bool allowRename )
ufkt.integral.lineWidth = n.attribute( "integral-width" ).toDouble() * lengthScaler;
ufkt.integral_use_precision = n.attribute( "integral-use-precision" ).toInt();
ufkt.integral_precision = n.attribute( "integral-precision" ).toInt();
ufkt.startx.updateExpression( n.attribute( "integral-startx" ) );
ufkt.starty.updateExpression( n.attribute( "integral-starty" ) );
ufkt.setIntegralStart( n.attribute( "integral-startx" ), n.attribute( "integral-starty" ) );
QDomElement minElement = n.namedItem( "arg-min" ).toElement();
QString expression = minElement.text();
......
......@@ -81,8 +81,7 @@ k_dcop:
virtual bool setFunctionMaxValue(const QString &max, uint id) = 0;
virtual QString functionStartXValue(uint id) = 0;
virtual QString functionStartYValue(uint id) = 0;
virtual bool setFunctionStartXValue(const QString &x, uint id) = 0;
virtual bool setFunctionStartYValue(const QString &y, uint id) = 0;
virtual bool setFunctionStartValue(const QString &x, const QString &y, uint id) = 0;
};
......
......@@ -139,12 +139,12 @@ bool XParser::getext( Function *item, const QString fstr )
return true;
}
double XParser::a1fkt( Equation *u_item, double x, double h )
double XParser::derivative1( Equation *u_item, double x, double h )
{
return ( fkt(u_item, x + h ) - fkt( u_item, x ) ) / h;
}
double XParser::a2fkt( Equation *u_item, double x, double h )
double XParser::derivative2( Equation *u_item, double x, double h )
{
return ( fkt( u_item, x + h + h ) - 2 * fkt( u_item, x + h ) + fkt( u_item, x ) ) / h / h;
}
......@@ -247,12 +247,34 @@ void XParser::fixFunctionName( QString &str, Equation::Type const type, int cons
}
}
double XParser::euler_method(const double x, Equation * eq )
double XParser::integral( Equation * eq, double x, double h )
{
double const y = eq->oldy + ((x-eq->oldx) * eq->oldyprim);
eq->oldy = y;
eq->oldx = x;
eq->oldyprim = fkt( eq, x ); //yprim;
h = qAbs(h);
assert( h > 0 ); // in case anyone tries to pass us a zero h
// the difference between h and dx is that h is only used as a hint for the
// stepwidth; dx is made similar to h in size, yet tiles the gap between x
// and the previous x perfectly
// we use the 2nd degree Newton-Cotes formula (Simpson's rule)
double a0 = eq->lastIntegralPoint.x();
double y = eq->lastIntegralPoint.y();
int intervals = qRound( qAbs(x-a0)/h );
double dx = (x-a0) / intervals;
for ( int i = 0; i < intervals; ++i )
{
double a = a0 + (dx*i);
double m = a + (dx/2);
double b = a + dx;
y += (dx/6.0)*( fkt( eq, a ) + (4.0 * fkt( eq, m )) + fkt( eq, b ) );
}
eq->lastIntegralPoint = QPointF( x, y );
return y;
}
......@@ -518,36 +540,28 @@ bool XParser::setFunctionMaxValue(const QString &max, uint id)
return true;
}
QString XParser::functionStartXValue(uint id)
{
if ( !m_ufkt.contains( id ) )
return 0;
return m_ufkt[id]->startx.expression();
}
bool XParser::setFunctionStartXValue(const QString &x, uint id)
bool XParser::setFunctionStartValue(const QString &x, const QString &y, uint id)
{
if ( !m_ufkt.contains( id ) )
return false;
m_ufkt[id]->startx.expression() = x;