forum.alglib.net
http://forum.alglib.net/

spare matrix
http://forum.alglib.net/viewtopic.php?f=2&t=3855
Page 1 of 1

Author:  Emmacalculus007 [ Thu May 24, 2018 2:10 pm ]
Post subject:  spare matrix

hello admin, please help check my C++ code. it complied but didnt return 0, so n output.
The codes are below;
//HOWARTH TRANSFORMATION OF BLASIUS BOUNDARY LAYER PRROBLEM
// USING 100STEPS i.e t=0.05,
//ASSUMING INFINITY=5

#include<iostream>
#include "stdafx.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "solvers.h"
#include<cmath>
#include<iomanip>

using namespace std;
using namespace alglib;


int main(int argc, char **argv)
{

int m,n,i,j,k,M;// sum={0};
double A[300][1], B[300][300], C[300][1], D[300][1], E[300][1], F[300][1], s[300]={0}, f[300], h[300], T[300] ,t,a,b,w,z,y,x, p;// l[300][300]={0},u[300][300]={0},z[300]={0};


cout<<"Aim: To solve Howarth Transformation Boundary Layer problem using nonlinear solver (Newton Raphson method) for: \n"<<" f'''+ f/(h^(-1/3)) f''= 0 \n"<<" h''+ (0.7fh'/h^(-1/3))- 0.28M^(2)*f''^(2)=0 \n"<<endl;

cout<<"Enter number of unknown variables "<<endl;
cin>>m;
cout<<"Enter number of non-linear equations "<<endl;
cin>>n;
cout<<"Enter initial guess value for f "<<endl;
cin>>a;
cout<<"Enter initial guess value of for h "<<endl;
cin>>b;

for(i=1;i<=99;i++)
{
f[i]= a;
h[i]= b;
}



/* 198 non-linear functions with 198 unknowns Fn(f1,f2,f3,f4,h1,h2,h3,h4) */

A[0][0]= -2*f[1] - 2*f[2] + f[3] + (0.1*f[1]*f[2])/pow(h[1],-1/3) - ((0.2*pow(f[1],2))/pow(h[1],-1/3));
A[1][0]= 2*f[1] - 2*f[3] + f[4] + (0.1*f[1]*f[2])/pow(h[2],-1/3) + (0.1*f[2]*f[3])/pow(h[2],-1/3) - (0.2*pow(f[2],2))/pow(h[2],-1/3);

for(i=3;i<=97;i++)
{
A[i-1][0]= f[i+2]-2*f[i+1]+2*f[i-1]-2*f[i-2] + 0.1*f[i]*f[i+1]/pow(h[i],-1/3) + 0.1*f[i]*f[i-1]/pow(h[i],-1/3) - (0.2*pow(f[i],2))/pow(h[i],-1/3);
}

A[97][0]= 0.1-f[99]+ 2*f[97]-2*f[96] + (0.1*f[98]*f[99])/pow(h[98],-1/3) +(0.1*f[98]*f[97])/pow(h[98],-1/3) - ((0.2*pow(f[98],2))/pow(h[98],-1/3));
A[98][0]= -f[99]+ 2*f[98]-2*f[97] + (0.005*f[99])/pow(h[99],-1/3) + (0.1*f[99]*f[98])/pow(h[99],-1/3) - (0.1*pow(f[99],2))/pow(h[99],-1/3) + 0.05;
A[99][0]= pow(M,2)*(-2.24*pow(f[1],2) - 0.56*pow(f[2],2) + 2.24*f[1]*f[2] )- 0.01*h[1] + 0.005*h[2] + (0.0000875*f[1]*h[2])/pow(h[1],-1/3) - (0.000175*f[1])/pow(h[1],-1/3) + 0.01;

for(i=2; i<=98; i++)
{
A[98+i][0]= pow(M,2)*(-2.24*pow(f[i],2) - 0.56*pow(f[i-1],2) - 0.56*pow(f[i+1],2) + 2.24*f[i]*f[i-1] - 1.12*f[i-1]*f[i+1] + 2.24*f[i]*f[i+1] )- 0.01*h[i] + 0.005*h[i-1] + 0.005*h[i+1] + (0.0000875*f[i]*h[i+1])/pow(h[i],-1/3) - (0.0000875*f[i]*h[i-1])/pow(h[i],-1/3);
}

A[197][0]= pow(M,2)*(-0.0014 -0.0056*f[98] + 0.0056*f[99] - 0.56*pow(f[98],2) - 0.56*pow(f[99],2) + 1.12*f[98]*f[99]) + 0.005*h[98] - 0.01*h[99] + (0.0000875*f[99])/pow(h[99],-1/3) - (0.0000875*f[99]*h[98])/pow(h[99],-1/3) + 0.005 ;


/* Jacobian Matrix elements of Matrix B */

B[0][0]= -2+(0.1*f[2])/pow(h[1],-1/3)-(0.4*f[1])/pow(h[1],-1/3); B[0][1]= -2+(0.1*f[1])/pow(h[1],-1/3); B[0][2]= 1; B[0][3]= 0; B[0][99]= (0.1*f[1]*f[2]*pow(h[1],-2/3)/3)-(0.2*pow(f[1],2)*pow(h[1],-2/3)/3); B[0][100]= 0; B[0][101]= 0; B[0][102]= 0;

B[1][0]= 2+(0.1*f[2])/pow(h[2],-1/3); B[1][1]= (0.1*f[1])/pow(h[2],-1/3)+ (0.1*f[3])/pow(h[2],-1/3)-(0.4*f[2])/pow(h[2],-1/3); B[1][2]= -2+ (0.1*f[2])/pow(h[2],-1/3) ; B[1][3]= 1; B[1][99]= 0; B[1][100]= (0.1*f[1]*f[2]*pow(h[2],-2/3)/3) + (0.1*f[2]*f[3]*pow(h[2],-2/3)/3) - (0.2*pow(f[2],2)*pow(h[2],-2/3)/3); B[1][101]= 0; B[1][102]= 0;

for(i=3; i<=97; i++)
{
B[i-1][i-1]= 0.1*f[i+1]/pow(h[i],-1/3) - 0.4*f[i]/pow(h[i],-1/3) + 0.1*f[i-1]/pow(h[i],-1/3);
B[i-1][i-2]= 2 + 0.1*f[i]/pow(h[i],-1/3);
B[i-1][i+1] = 1 ;
B[i-1][i]= -2 + 0.1*f[i]/pow(h[i],-1/3);
B[i-1][i-3] = -2 ;
B[i-1][98+i]= (0.1*f[i]*f[i+1]*(h[i],-2/3)/3) - (0.2*(f[i],2)*(h[i],-2/3)/3) + (0.1*f[i]*f[i-1]*(h[i],-2/3)/3);
}

B[97][95]= -2; B[97][96]= 2+(0.1*f[98])/pow(h[98],-1/3); B[97][97]= (0.1*f[99])/pow(h[98],-1/3) + (0.1*f[97])/pow(h[98],-1/3)- (0.4*f[98])/pow(h[98],-1/3); B[97][98]= -1 + (0.1*f[98])/pow(h[98],-1/3); B[97][196]= (0.1*f[98]*f[99]*pow(h[98],-2/3)/3) + (0.1*f[98]*f[97]*pow(h[98],-2/3)/3)-(0.2*pow(f[98],2)*pow(h[98],-2/3)/3);

B[98][96]= -2; B[98][97]= 2+(0.1*f[99])/pow(h[99],-1/3); B[98][98]= -1+ (0.005)/pow(h[99],-1/3)+ (0.1*f[98])/pow(h[99],-1/3)-(0.2*f[99])/pow(h[99],-1/3); B[98][197] = (0.1*f[99]*f[98]*pow(h[99],-2/3)/3) + (0.01*f[99]*pow(h[99],-2/3)/3)-(0.1*pow(f[99],2)*pow(h[99],-2/3)/3);

B[99][0]= pow(M,2)*(-4.48*f[1] + 2.24*f[2]) + (0.0000875*h[2])/pow(h[1],-1/3) - (0.000175)/pow(h[1],-1/3) ; B[99][1]= pow(M,2)*(2.24*f[1] - 1.12*f[2]); B[99][2]= 0; B[99][3]= 0; B[99][99]= -0.01 + (0.0000875*f[1]*h[2]*pow(h[1],-2/3)/3) - ( 0.000175*f[1]*pow(h[1],-2/3)/3); B[99][100]= 0.005 + (0.0000875*f[1])/pow(h[1],-1/3); B[99][101]= 0; B[99][102]= 0;

for(i=2; i<=98; i++)

{
B[98+i][i-1]= pow(M,2)*(2.24*f[i+1] - 4.48*f[i] + 2.24*f[i-1] ) + (0.0000875*h[i+1])/pow(h[i],-1/3) - (0.0000875*h[i-1])/pow(h[i],-1/3) ;
B[98+i][i-2]= pow(M,2)*(-1.12*f[i+1] + 2.24*f[i] - 1.12*f[i-1]);
B[98+i][i]= pow(M,2)*(-1.12*f[i+1] + 2.24*f[i] - 1.12*f[i-1] );
B[98+i][98+i]= -0.01 + (0.0000875*f[i]*h[i+1]*pow(h[i],-2/3)/3 ) - (0.0000875*f[i]*h[i-1]*pow(h[i],-2/3)/3 );
B[98+i][97+i]= 0.005 - (0.0000875*f[i])/pow(h[i],-1/3) ;
B[98+i][99+i]= 0.005 + (0.0000875*f[i])/pow(h[i],-1/3) ;
}

B[197][97]= pow(M,2)*(-0.0056 - 1.12*f[98] + 1.12*f[99] ); B[197][98]= pow(M,2)*(0.0056 - 1.12*f[99] + 1.12*f[98] ) + (0.0000875/pow(h[99],-1/3) - (0.0000875*h[98])/pow(h[99],-1/3)); B[197][196]= 0.005 - (0.0000875*f[99])/pow(h[99],-1/3) ; B[197][197]= -0.01 + (0.0000875*f[99]*pow(h[99],-2/3)/3) - ( 0.0000875*f[99]*h[98]*pow(h[99],-2/3)/3);


//for(i=0;i<198;i++)
// for(j=0;j<198;j++)
// {
// ae_int_t w = 198;
// sparsematrix B[i][j];
// real_1d_array z= "A[i][j]";
// real_1d_array s[i];
// sparsesolverreport rep;
// bool isuppertriangle = false;
// sparsesolvesks(B[i][j], w, isuppertriangle, z, rep, s[i]);
//}
// printf("%s\n", s[i]);
//



for(i=0;i<198;i++)
for(j=0;j<198;j++)
{

real_2d_array y= "B[i][j]";
ae_int_t n = 198;
bool isupper = false;
real_2d_array z = "A[i][j]";
ae_int_t m = 198;
ae_int_t info;
densesolverreport rep;
real_2d_array s[i];
smp_spdmatrixsolvem(y, n, isupper, z, m, info, rep, s[i]);
}
printf("%s\n", s[i]);
return 0;
}

Page 1 of 1 All times are UTC
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group
http://www.phpbb.com/