forum.alglib.net

ALGLIB forum
 It is currently Fri Sep 29, 2023 10:52 am

 All times are UTC

Forum rules

1. This forum can be used for discussion of both ALGLIB-related and general numerical analysis questions
2. This forum is English-only - postings in other languages will be removed.

 Page 1 of 1 [ 1 post ]
 Print view Previous topic | Next topic
Author Message
 Post subject: spare matrixPosted: Thu May 24, 2018 2:10 pm

Joined: Mon May 21, 2018 8:07 pm
Posts: 1
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;
}

Top

 Display posts from previous: All posts1 day7 days2 weeks1 month3 months6 months1 year Sort by AuthorPost timeSubject AscendingDescending
 Page 1 of 1 [ 1 post ]

 All times are UTC

Who is online

Users browsing this forum: No registered users and 1 guest

 You cannot post new topics in this forumYou cannot reply to topics in this forumYou cannot edit your posts in this forumYou cannot delete your posts in this forumYou cannot post attachments in this forum

Search for: