Trapezoid multiparameter function

[zaahm18]

Hi, I have an assignment to deliver until next Monday.

I wanted to use this function:
F(x,y,d) = 1 / 2 * (1 / pow((d * d + (x - y) * (x - y)), 3 / 2)


(C++ function:

Doub Func(Doub x, Doub y, Doub d) {
	return 1 / 2 * (1 / pow((d * d + (x - y) * (x - y)), 3 / 2));
}

to be used by the trapezoid integration method:

template<class T>
Doub trapezoidal(T& f, const Doub a, const Doub b, const int n) {
	if (n < 2) throw("n must be at least 2 for trapezoidal");
	Doub h = (b - a) / (n - 1);
	Doub sum = 0.5 * (f(a) + f(b));
	Doub x = a + h;
	for (int i = 2; i < n; i++) {
		sum += f(x);
		x += h;
	}
	return sum * h;
}

I use also the functional library to define function types:

typedef function<Doub(Doub, Doub, Doub)> F;
typedef function<Doub(F&)> F2;
typedef function<Doub(F2&, const Doub, const Doub, const int)> M;

But somehow I can't implement the function Func into the trapezoid function. how do I do it?

In the case of Doub, replace with double, because I use numeric methods include files.


The reason why I cant use trapezoid(Func,a,b,n), because Func uses 3 parameters, and my trapezoid takes function of only one.

[lastchance]

As it stands your function is going to return 0 every time, because of integer division.

[zaahm18]

No. None of them are integers. They all are doubles.

My version of double is Doub, which uses the NR.3 version (numerical recipe).

[seeplus]

But in C++, 3 / 2 is 1 and 1 / 2 is 0 as integer division is used. You need:

Doub Func(Doub x, Doub y, Doub d) {
	return 0.5 * (1.0 / pow((d * d + (x - y) * (x - y)), 1.5));
}

[seeplus]

Your function f() is defined in 2 variables x and y (d is a constant?)

In maths, how would you use the trapezoid rule to integrate a function of 2 variables over a stated range? Once you know how to do this in maths using pen/paper/calculator then you can translate that algorithm into a program design and then into code.

Sorry but my calculus only extends to functions of 1 variable.

[zaahm18]

Ok, I've modified the integers.

and yes, d is a constant, and one of the variables of the function (x or y) will remain the constant:

here is the pseudocode for the actual integration for the Func:

u(x) = 1.0+(1-0.5)*Int(f(x,y,d),-w/2,w/2,dy)

I want to use the trapezoid method to integrate a function, which has three parameters, 2 of which are constants (the example earlier I've shown is that, y and d are constants, while x is the variable)

This is the main function I want to run through:

	

int main()
{
F fs = Func; 
M method = trapezoidal<F2>;
Doub trapezoidal_k = 2;
vector<string> plate_names = { "Plate 1", "Plate 2" };
// the vector of given values N:
vector<Int> _N = { 3,6, 20,56,100 };
	
// This range based for loop runs entirely on the vector of given values.
// problem c) and d) are solved inside this for loop:
for (auto N : _N) {
// c:
// This for-loop initializes all values for both A1 and A2 matrices:
// And create the design matrix A out of it.
MatDoub A1(N+1, N+2), A2(N+1,N+2);
Doub x, y, Aj1, Aj2;
F2 F_1, F_2;
int n = 1;
    for (int i = 0; i <= N; i++)
	{
	         A1[i][0] = epsilon1 * sigma * pow(T1, 4);
	         A2[i][0] = epsilon2 * sigma * pow(T2, 4);
	         for (int j = 0; j <= N; j++)
		{
			Aj1 = fs(3,4,d);
			Aj2 = fs(3,4,d);
			Aj1 = method(F_1, -(1 / 2 * w), (1 / 2 * w), n);
			Aj2 = method(F_2, -(1 / 2 * w), (1 / 2 * w), n);

			A1[i][j + 1] = (1-epsilon1)*Aj1;
			A2[i][j + 1] = (1-epsilon2)*Aj2;
		}
	}
}

[zaahm18]

Note that I've used typedef to cut down the parameters of Func from 3 to 1.

[lastchance]

OK, so A1[i][0] and A2[i][0] contain radiative heat transfer (Stefan-Bolzmann law) from two surfaces at uninitialised absolute temperatures T1 and T2, ... using emissivities epsilon1 and epsilon2 and a Stefan-Bolzmann constant sigma that also isn't defined in the code you've shown ...

Then you set Aj1 and Aj2 to ... something or other ... which you immediately overwrite with some other stuff ...

... which still isn't right, because you ignored my first post and still haven't realised that 1 / 2 will give 0 in c++ (and many other languages) because of integer division ... so you won't be integrating from -w/2 to w/2 but from 0 to 0 ...


It would be quicker if you just stated your assignment.

[zaahm18]

Sorry, i was in a rush:

but anyway, here are the constants:

Int T1 = 1000, T2 = 500;
Doub epsilon1 = 0.8, epsilon2 = 0.6, sigma = 1.712e-9, d = 1.0, w = 1.0;

The assignment is divided in 4 section:

a), b), c), d)

a)
Here we have to approximate the both integrals Int(F(x,y,d)v(y),-w/2,w/2,dy) and Int(F(x,y,d)u(x),-w/2,w/2,dx) using trapexoid method, and we should write it on paper. I have done this part. (index goes from 0 to N, not N-1, there are N+1 indexes in this case).

these are the equations u(x) = epsilon1*sigma*T1^4+(1-epsilon1)*Int(F(x,y,d)v(y),-w/2,w/2,dy).

same with v(y) but with 2 and opposite integral.

b)
Here we use the derived part from a) and write it into the discrete problem, so it becomes the system of linear equation.

c)
this is the part I am stuck in Write the program that sets up and solves that linear equation from b) for the given value N (this case is 4). Although I know how to solve it, but that is only possible if I can somehow integrate Func using the trapezoid method.

[lastchance]

What do you mean by "the system of linear equation"?

The most obvious interpretation is that u_i is written as a weighted sum of the v_i (and vice versa), with the weighted sum replacing the integral. In which case you never actually have to work out the integral by the trapezium rule.

Your code bears no resemblance to the problem.

[zaahm18]

That because the first two problems aren't supposed to solve with code.

and the second problem, you will get a system of linear equation:

eps1*sig*T1^4 + (1-eps1)*h*(1/2 F(x,y0,d)*v0+F(x,y1,d)v1+...+1/2F(x,yN,d)vN)

from that equation, you create a design matrix (but that is not the issue here)

And the c) problem is where you solve the linear equation from b), but that part I am not having trouble with. The part of integrating F(x,y,d), by using trapezoid method is where I am stuck.

[lastchance]

Your linear equation above IS the algebraic statement of the trapezium rule.
It relates each u_i to v0, v1, ...., vn (and vice versa)
This will be 2(N+1) equations in total, which you can solve by any matrix-inversion technique, NOT by integration by trapezium rule.

None of your stated problems actually require you to carry out a trapezium-rule integration in code.

[zaahm18]

look, this is my trapezoid method:

template<class T>
Doub trapezoidal(T& f, const Doub a, const Doub b, const int n) {
	if (n < 2) throw("n must be at least 2 for trapezoidal");
	Doub h = (b - a) / (n - 1);
	Doub sum = 0.5 * (f(a) + f(b));
	Doub x = a + h;
	for (int i = 2; i < n; i++) {
		sum += f(x);
		x += h;
	}
	return sum * h;
}

and this is my function I have given:

Doub Func(Doub x, Doub y, Doub d) {
	return 1 / 2 * (1 / pow((d * d + (x - y) * (x - y)), 3 / 2));
}

I almost forgot to tell you that; the function Func, which takes 3 parameters, will only have 1 parameter with unknown because d is a constant, and one of the variables will have a value from the definite integral, only one variable will be unknown, and that will be the function variable:

suppose you have Int(F(x,y,d),1,3,dx);

the variable x will go from 1 to 3 after finding its antiderivative, and the d is the constant.

We only get a function that takes 1 parameter instead of 3 this time.

[lastchance]

Please upload your COMPLETE ORIGINAL assignment - NOT YOUR RENDITION of it. This is rapidly becoming an XY problem: https://en.wikipedia.org/wiki/XY_problem

Your trapezoidal integration function wasn't given to you - it's just taken from "Numerical Recipes". It is daft here to use "Doub".

I repeat:
- there is nothing in what you have actually stated of your problem thus far (you gave your rendition of parts (a)-(c), but not (d)) which suggests that you actually need to do numerical integration in code;
- (as @Seeplus has also told you) if you write 1 / 2 in c++ you will get 0 and if you write 3 / 2 you will get 1, both because of the truncation incurred by integer division.

Please don't keep repeating your 3-parameter/1-parameter dirge - just state the exact, original problem.

[zaahm18]

The assignment is in pdf file

file:///C:/Users/zaina/OneDrive/Skoleopgaver/Robotteknologi/Semester_6/Nummerisk%20metoder/Assignment/Assigment_2_plate_radiation/plate_radiation.pdf

(you can't open it)

Two infinitely long parallel plates of width w are placed with distanced.
The plates radiate and reflect diffusely and have uniform temperatures T1 and T2, and
emmisivities ε1 and ε2, respectively. Let u(x) denote the radiation per unit time per unit
area from plate 1 at the distance x from the mid axis of the plate, and let v(y) be the
corresponding function for plate 2. Then u and v are solutions to the coupled integral
equations:

u(x) = ε1σT₁⁴+ (1 − ε1)*Int(F(x,y,d)v(y),-1/2w,1/2w,dy)

v(y) = ε2σT₂⁴+ (1 − ε2)*Int(F(x,y,d)u(x),-1/2w,1/2w,dx)


where σ is the Stefan-Bolzmann constant, and F is given by

F(x,y,d) = 1/2*1/(d²+(x+y)²)^3/2

The integrals Int(F(x,y,d)v(y),-1/2w,1/2w,dy) and Int(F(x,y,d)u(x),-1/2w,1/2w,dx)

describe the incoming radiation per unit time per unit area from the other plate. The
total energy (per unit time and unit length).

Data: T₁ = 1000, T₂ = 500, ε₁ = 0.80, ε₂ = 0.60, σ = 1.712 · 10⁻⁹
, d = 1.00, w = 1.00.


a) Use the trapezoidal method to write an approximation of (4) and (5). Hint: Introduce (u_i)_i=0^N and (v_i)_i=0^N describing the values of u and v in the points (x_i)^N_i=0 and (y_i)^N_i=0, respectively.

b) Using the approximation from a), write the discrete problem corresponding to (1)
and (2) as a system of linear equations in the variables (u_i)_i=0^N and (v_i)_i=0^N

c) Write a program that sets up and solves the system of linear equations from b) for
a given value of N. N = 4:

State the found values of u(x) for x = −1/2, −1/4, 0,1/4,1/2 and v(y) for y = −
1/2, −1/4, 0,1/4,1/2.
Hint: Each u(x) and v(y), you’re asked to state, correspond to one of the variables (u_i)^N_i=0, (v_i)^N_i=0.

[lastchance]

#include <iostream>
#include <vector>
#include <cmath>
using namespace std;

using vec    = vector< double >;
using matrix = vector< vec >;


//======= Function prototypes ========================================

bool solve( matrix A, vec B, vec &X );
double F( double x, double y, double d );        // The infamous function

//====================================================================

int main()
{
   const double epsilon1 = 0.8, epsilon2 = 0.6;  // Emissivities
   const double sigma = 1.712e-9;                // Stefan Boltzmann constant
   const double T1 = 1000.0, T2 = 500.0;         // Absolute temperatures of the plates
   const double d = 1.0;                         // Separation of the plates
   const double w = 1.0;                         // Width of the plates

   const int N = 4;                              // Number of intervals
   int NN = 2 * ( N + 1 );                       // Number of degrees of freedom (values of u and v)

   // Set up coordinates
   vec x(1+N);                
   double dx = w / N;
   for ( int i = 0; i <= N; i++ ) x[i] = -0.5 * w + i * dx;
   vec y = x;
   double dy = dx;

   // Set up linear system as (there are other possibilities)
   //  M {u} = b
   //    {v}
   matrix M( NN, vec(NN,0.0) );
   vec B(NN);
   vec XX(NN);
   double heat1 = epsilon1 * sigma * T1 * T1 * T1 * T1;
   double heat2 = epsilon2 * sigma * T2 * T2 * T2 * T2;
   for ( int i = 0; i <= N; i++ )                          // u equation (first N+1 rows of M; x = x[i])
   {
      B[i] = heat1;
      M[i][i] = 1.0;                                       // u[i] = ...
      M[i][N +1] = -( 1.0 - epsilon1 ) * dy * 0.5 * F( x[i], y[0], d );        // "Integral" transferred to LHS
      M[i][NN-1] = -( 1.0 - epsilon1 ) * dy * 0.5 * F( x[i], y[N], d );
      for ( int j = 1; j < N; j++ ) M[i][N+1+j] = -( 1.0 - epsilon1 ) * dy * F( x[i], y[j], d );
   }
   for ( int i = 0; i <= N; i++ )                          // v equation (last N+1 rows of M; y = y[i])
   {
      B[N+1+i] = heat2;
      M[N+1+i][N+1+i] = 1.0;                               // v[i] = ....
      M[N+1+i][0] = -( 1.0 - epsilon2 ) * dx * 0.5 * F( x[0], y[i], d );
      M[N+1+i][N] = -( 1.0 - epsilon2 ) * dx * 0.5 * F( x[N], y[i], d );
      for ( int j = 1; j < N; j++ ) M[N+1+i][j] = -( 1.0 - epsilon2 ) * dy * F( x[j], y[i], d );
   }

   // Solve
   cout << "Solving ... " << boolalpha << solve( M, B, XX ) << '\n';
   vec U( XX.begin(), XX.begin()+1+N );
   vec V( XX.begin()+1+N, XX.end()   );

   cout << "\nU values: ";   for ( double e : U ) cout << e << "  ";
   cout << "\nV values: ";   for ( double e : V ) cout << e << "  ";
}

//====================================================================


double F( double x, double y, double d )
{
   double rsq = d * d + ( x - y ) * ( x - y );
   return 0.5 / rsq / sqrt( rsq );
}


//====================================================================


bool solve( matrix A, vec B, vec &X )
//--------------------------------------
// Solve AX = B by Gaussian elimination
//--------------------------------------
{
   const double SMALL = 1.0e-50;

   int n = A.size();

   // Row operations for i = 0, ,,,, n - 2 (n-1 not needed)
   for ( int i = 0; i < n - 1; i++ )
   {
      // Pivot: find row r below with largest element in column i
      int r = i;
      double maxA = abs( A[i][i] );
      for ( int k = i + 1; k < n; k++ )
      {
         double val = abs( A[k][i] );
         if ( val > maxA )
         {
            r = k;
            maxA = val;
         }
      }
      if ( r != i )
      {
         for ( int j = i; j < n; j++ ) swap( A[i][j], A[r][j] );
         swap( B[i], B[r] );
      }

      // Row operations to make upper-triangular
      double pivot = A[i][i];
      if ( abs( pivot ) < SMALL ) return false;            // Singular matrix

      for ( int r = i + 1; r < n; r++ )                    // On lower rows
      {
         double multiple = A[r][i] / pivot;                // Multiple of row i to clear element in ith column
         for ( int j = i; j < n; j++ ) A[r][j] -= multiple * A[i][j];
         B[r] -= multiple * B[i];
      }
   }
   if ( abs( A[n-1][n-1] ) < SMALL ) return false;         // Singular matrix

   // Back-substitute
   X = B;
   for ( int i = n - 1; i >= 0; i-- )
   {
      for ( int j = i + 1; j < n; j++ )  X[i] -= A[i][j] * X[j];
      X[i] /= A[i][i];
   }

   return true;
}


//==================================================================== 

Solving ... true

U values: 1390.15  1394.08  1395.6  1394.08  1390.15  
V values: 260.505  297.021  311.005  297.021  260.505  

[zaahm18]

Tnx for the suggestion ;).

[zaahm18]

I am maybe stuck in whole another problem

:Use the results to find Q1 and Q2.

where:

Q1 = Int(u(x)-I1(x),-1/2w,1/2w,dx).

and

Q2 = Int(v(y)-I2(y),-1/2w,1/2w,dy).

my result are:

Q1: 1361.12 and Q2: 138.925
Q1: 1376.61 and Q2: 212.555
Q1: 1385.14 and Q2: 252.76
Q1: 1389.59 and Q2: 273.705
Q1: 1391.85 and Q2: 284.387
Q1: 1393 and Q2: 289.781
Q1: 1393.57 and Q2: 292.491

instead of I use N = 4, I use N = {4,8,16,32,64,128,256}, so that's why I get 7 different results.

this is my code:

		
// We use result to find Q1 and Q2.
Doub Q1 = 0, Q2 = 0, I1 = 0, I2 = 0;

// Define u(x) and v(y) (u and v):
VecDoub u(N + 1), v(N + 1);
for (int i = 0; i <= N; i++) {
	u[i] = X[i];
	v[i] = X[N + 1 + i];
}

// This for-loop runs the integration of (ux or vy - I1 or I2):
for (int i = 0; i <= N; i++)
{
	if (i == 0)
	{
		I1 = v[i] * dy * 0.5 * F(x[i], y[0], d);
		I2 = u[i] * dx * 0.5 * F(x[0], y[i], d);
		Q1 += dx * 0.5 * (u[i] - I1);
		Q2 += dy * 0.5 * (v[i] - I2);
	}
	else if (i == N)
	{
	I1 = v[i] * dy * 0.5 * F(x[i], y[N], d);
	I2 = u[i] * dx * 0.5 * F(x[N], y[i], d);
	Q1 += dx * 0.5 * (u[i] - I1);
	Q2 += dy * 0.5 * (v[i] - I2);
	}
	else
	{
	I1 = v[i] * dy * F(x[i], y[i], d);
	I2 = u[i] * dx * F(x[i], y[i], d);
	Q1 += dx * (u[i] - I1);
	Q2 += dy * (v[i] - I2);
	}
}

cout << "The values of Q1 is: " << Q1 << " and Q2: " << Q2 << endl;
spaces();
line();

[lastchance]

You need to give a clearer statement of the problem. NOT in code.

[zaahm18]

its the task C:

c) Write a program that sets up and solves the system of linear equations from b) for a given value of N. Uses the results to find Q1 and Q2.
For each value of N = {4, 8, 16, 32, 64, 128, 256}:


State the found values of u(x) for x = {−1/2, −1/4, 0,1/4,1/2} and v(y) for y = {−1/2, −1/4, 0,1/4,1/2}.
Hint: Each u(x) and v(y), you’re asked to state, correspond to one of the variables (u_i)^N_i=0, (v_i)^N_i=0.

State Q₁ and Q₂.

I have already found u(x) and v(y).

but now I have to find I₁(x) and I₂(y).

The integrals which describe I₁(x) and I₂(y) are:

I₁(x) = Int(F(x,y,d)v(y),-1/2w,1/2w,dy)

I₂(y) = Int(F(x,y,d)u(x),-1/2w,1/2w,dx)

I believe we can fetch that value somewhere...

and the Qs:

Q₁ = Int(u(x)-I₁(x),-1/2w,1/2w,dx).

and

Q₂ = Int(v(y)-I₂(y),-1/2w,1/2w,dy).

my colleague told me that Q2 should have negative values, but I only get the positives, I've tried to debug my code, but I couldn't find the fault.