how to plot curves with complex numbers?

[rambo83]

Hello,

I have ploted the same functions in MATLAB and then with QwtPlot and the results are different. Could you tell me why and how to deal properly with complex numbers in QwtPlot?

Thank you

best regards,

Vitali

here is the definition of functions:

complex<double> benchfunc_2_constr_1(complex<double> x){
						return (5.0 + sqrt(100.0 - pow(x - 5.0 ,2)));			// (-pow(x1-5,2) - pow(x2-5,2) +100) <=0
						}												// resp. y >= 5 + sqrt(100 - pow(x - 5,2))
 
				complex<double> benchfunc_2_constr_2(complex<double> x){
						return (5.0 + sqrt(82.81 - pow(x - 6.0 ,2)));			// (-pow(x1-6,2) - pow(x2-5,2) - 82.81) <=0
						}												// resp. y <= 5 + sqrt(82,81 - pow(x - 6,2))

complex<double> benchfunc_2_constr_1(complex<double> x){
						return (5.0 + sqrt(100.0 - pow(x - 5.0 ,2)));			// (-pow(x1-5,2) - pow(x2-5,2) +100) <=0
						}												// resp. y >= 5 + sqrt(100 - pow(x - 5,2))

				complex<double> benchfunc_2_constr_2(complex<double> x){
						return (5.0 + sqrt(82.81 - pow(x - 6.0 ,2)));			// (-pow(x1-6,2) - pow(x2-5,2) - 82.81) <=0
						}												// resp. y <= 5 + sqrt(82,81 - pow(x - 6,2))

To copy to clipboard, switch view to plain text mode 

here is the class which generates data for ploting the curves:

class SimpleData: public QwtData
{
 
public:
    SimpleData( complex<double> (*b)(complex<double>), double rang1[2], double rang2[2], double step)
    {
    	constraintFunc = b;
    	_step = step;
    	rangeX[0]= rang1[0];	rangeX[1] = rang1[1];
    	rangeY[0]= rang2[0];	rangeY[1] = rang2[1];
    	range = (rangeX[1] - rangeX[0])/step;
    	if(range<0){
    		range = range *(-1);
    	}
    	d_size = size_t(range);
    }
 
    virtual QwtData *copy() const
    {
        return new SimpleData(*this);
    }
 
    virtual size_t size() const
    {
        return d_size;
    }
 
    virtual double x(size_t i) const
    {
    	// here the array with x values will be produced
    	return rangeX[0] + i*((rangeX[1]- rangeX[0])/range);
    }
 
    virtual double y(size_t i) const
    {
    	complex<double> result = (*constraintFunc)(x(i));	// compute the y value
 
    	if(real(result) > rangeY[1]){	// if exceeds the max bound
    		return rangeY[1];
    	}
    	else if(real(result) < rangeY[0]){ // if exceeds the min bound
    		return rangeY[0];
    	}
    	else return real(result);	// return only the real value of complex number
    }
private:
    size_t d_size;
    double range;
    double _step;
    complex<double> (*constraintFunc)(complex<double>);	// function pointer for constraint function
	double rangeX[2];
	double rangeY[2];
};

class SimpleData: public QwtData
{

public:
    SimpleData( complex<double> (*b)(complex<double>), double rang1[2], double rang2[2], double step)
    {
    	constraintFunc = b;
    	_step = step;
    	rangeX[0]= rang1[0];	rangeX[1] = rang1[1];
    	rangeY[0]= rang2[0];	rangeY[1] = rang2[1];
    	range = (rangeX[1] - rangeX[0])/step;
    	if(range<0){
    		range = range *(-1);
    	}
    	d_size = size_t(range);
    }

    virtual QwtData *copy() const
    {
        return new SimpleData(*this);
    }

    virtual size_t size() const
    {
        return d_size;
    }

    virtual double x(size_t i) const
    {
    	// here the array with x values will be produced
    	return rangeX[0] + i*((rangeX[1]- rangeX[0])/range);
    }

    virtual double y(size_t i) const
    {
    	complex<double> result = (*constraintFunc)(x(i));	// compute the y value

    	if(real(result) > rangeY[1]){	// if exceeds the max bound
    		return rangeY[1];
    	}
    	else if(real(result) < rangeY[0]){ // if exceeds the min bound
    		return rangeY[0];
    	}
    	else return real(result);	// return only the real value of complex number
    }
private:
    size_t d_size;
    double range;
    double _step;
    complex<double> (*constraintFunc)(complex<double>);	// function pointer for constraint function
	double rangeX[2];
	double rangeY[2];
};

To copy to clipboard, switch view to plain text mode 

somewhere in code the data for the curve is set with:

curve1->setData(SimpleData((*ii).second, rangex1, rangex2, 0.05));

curve1->setData(SimpleData((*ii).second, rangex1, rangex2, 0.05));

To copy to clipboard, switch view to plain text mode 

and here the MATLAB code for comparison:

[x1, x2] = meshgrid(13:1:100, 0:1:100);
f = (x1 -10).^3 + (x2 - 20).^3;    % problem function
g1 = @(x) sqrt(100 - (x-5).^2 ) + 5;            % constraint function 1
g2 = @(x) sqrt(82.81 -(x-6).^2 ) +5;            % constraint function 2
 
t=13:0.1:100;
y1 = g1(t);
y2 = g2(t);
 
 
figure
subplot(1,2,1)
surf(x1,x2,f);
subplot(1,2,2)
contour(x1,x2,f);
hold on 
plot(t, g1(t), 'r');
plot(t, g2(t), 'y');
hold off

[x1, x2] = meshgrid(13:1:100, 0:1:100);
f = (x1 -10).^3 + (x2 - 20).^3;    % problem function
g1 = @(x) sqrt(100 - (x-5).^2 ) + 5;            % constraint function 1
g2 = @(x) sqrt(82.81 -(x-6).^2 ) +5;            % constraint function 2

t=13:0.1:100;
y1 = g1(t);
y2 = g2(t);

    
figure
subplot(1,2,1)
surf(x1,x2,f);
subplot(1,2,2)
contour(x1,x2,f);
hold on 
plot(t, g1(t), 'r');
plot(t, g2(t), 'y');
hold off

To copy to clipboard, switch view to plain text mode 

[rambo83]

to define how the curves are "different" in MATLAB and QwtPlot: they look very similar if you not zoom in, but if you do, then you can see the difference, namely the curves intersect differently at x1 about 14...