#include <iostream>
#include <math.h>
#include "quanc8.h"
#include "splines.h"
using namespace std;
double x[100], f[100], xk[100], h, sx, lx, ux, A, B, FLAG, ERREST, RELERR, ABSERR, b[100],c[100],d[100], sR[100];
double xSpline, ss;
int k[100], NOFUN, i, n;
double FUN ( double t)
{
return exp(x*t)*sin(t);
}
void main()
{
cout<<"Insert lower limit of 'x'."<<endl; // начало промежутка для х, узлов
cin>>lx;
sx=lx;
cout<<"Insert upper limit of 'x'."<<endl; // конец промежутка для х, узлов
cin>>ux;
cout<<"Insert step of 'x'."<<endl; // шаг увеличения значения узлов
cin>>h;
cout<<"Insert lower limit of integration."<<endl; // нижний предел интегрирования для QUANC8
cin>>A;
cout<<"Insert upper limit of integration."<<endl; // верхний предел интегрирования для QUANC8
cin>>B;
cout<<"Insert absolut error."<<endl; // Абсолютная погрешность
cin>>ABSERR;
cout<<"Insert relative error."<<endl; // Относительная погрешность
cin>>RELERR;
cout<<"Value of 'X' | Value of F"<<endl;
//Вычисление и вывод значений функции в заданных узлах
for ( i=0; sx<=ux; i++ )
{
x=sx;
sx+=h;
quanc8( FUN, A, B, ABSERR, RELERR, &f, &ERREST, &NOFUN, &FLAG);
cout<<x<<"\t|\t"<<f<<endl;
}
cout<<"Insert step of 'x' for spline nodes."<<endl; // Шаг для узлов в которых необходимо вычислить сплайн-функцию
cin>>ss;
n=i;
i=0;
spline(n,x,f,b,c,d);
xSpline=0.1;
cout<<"Value of 'xSpline' | Value of function in area of 'xSpline'"<<endl;
// Вывод значений полученных в результате решения сплайна
for ( i=0; i<n; i++ )
{
sR=seval(n, &xSpline, x, f, b, c, d);
cout<<xSpline<<"\t|\t"<<sR<<endl;
xSpline+=ss;
}
}
Странная проблема с выводом double на консоль.
Здравствуйте. Вопрос действительно может показаться глупым.
При написании учебной программы для вычисления функции с помощью QUANC8 и SPLINE SEVAL столкнулся с небольшой проблемой функция seval вместо решения функции в окрестности заданной точки выдаёт непонятные значения:
Программа QUANC8:[INDENT]quanc8.h:
quanc8.cpp
[/INDENT]Код программ SPLINE, SEVAL:[INDENT]splines.h
splines.cpp
[/INDENT][/INDENT]
Среда разработки: Visual Studio 2008
Операционная система: Windows 7, а также проверенно на Windows XP
При написании учебной программы для вычисления функции с помощью QUANC8 и SPLINE SEVAL столкнулся с небольшой проблемой функция seval вместо решения функции в окрестности заданной точки выдаёт непонятные значения:
Цитата:
-1.#IND
Вроде проверил всё что может грешить не могу понять в чём дело, все подозрения падают на функцию seval.
Все кого спрашивал, говорят что проблема в индексации памяти, вопрос как решить эту проблему. Помогите, может быть кто-то знает как решить эту проблему или в чём ошибка? Очень на вас надеюсь
Прилагаю: [INDENT]Свой Код(Main.cpp):
Код:
Программа QUANC8:[INDENT]quanc8.h:
Код:
// quanc8.h
#ifndef _QUANC8
#define _QUANC8
void quanc8(double (*FUN)(double), double A, double B, double ABSERR,
double RELERR, double *RESULT, double *ERREST, int *NOFUN, double *FLAG);
#endif
#ifndef _QUANC8
#define _QUANC8
void quanc8(double (*FUN)(double), double A, double B, double ABSERR,
double RELERR, double *RESULT, double *ERREST, int *NOFUN, double *FLAG);
#endif
quanc8.cpp
Код:
// quanc8.cpp
#include "quanc8.h"
#include <math.h>
void quanc8(double (*FUN)(double), double A, double B, double ABSERR,
double RELERR, double *RESULT, double *ERREST, int *NOFUN, double *FLAG)
/*
Integrale di FUN(X) da A a B
Routine automatica di tipo adattivo basata
sulle formule di Newton-Cotes ad 8 intervalli.
Input
FUN nome della funzione integranda che viene fornita nella
function FUN(X) di tipo double con argomento double.
A limite inferiore di integrazione.
B limite superiore di integrazione.
RELERR errore relativo.
ABSERR errore assoluto.
Output
RESULT una approssimazione dell'integrale che dovrebbe soddisfare
la precisione data dal meno rigoroso dei due errori
ERREST stima della grandezza dell'errore commesso
NOFUN numero di volte in cui la funzione e' stata calcolata
FLAG indicatore di attendibilita'; se FLAG e' zero si e'
ottenuta la precisione richiesta.
se FLAG e' xxx.yyy allora xxx= il numero di intervalli in cui non e'
stata raggiunta la precisione richiesta, mentre 0.yyy ci da' una
indicazione sul punto nell'intorno del quale non si e' raggiunta la
precisione richiesta.
*/
{
double QRIGHT[32],F[17],X[17],FSAVE[9][31],XSAVE[9][31];
int LEVMIN,LEVMAX,LEVOUT,NOMAX,NOFIN,LEV,NIM,J,I;
double W0,W1,W2,W3,W4,COR11,AREA,X0,F0,STONE,STEP;
double QLEFT,QNOW,QDIFF,QPREV,TOLERR,ESTERR;
/* *** PASSO 1 *** INIZIALIZZAZIONE GENERALE COSTANTI */
LEVMIN=1;
LEVMAX=30;
LEVOUT=6;
NOMAX=5000;
NOFIN=NOMAX-(8*(LEVMAX-LEVOUT+pow(2,(double)(LEVOUT+1))));
/* ISTRUZIONI RIGUARDO AL CASO IN CUI NOFUN RAGGIUNGE IL VALORE NOFIN */
W0=3956.0/14175.0;
W1=23552.0/14175.0;
W2=-3712.0/14175.0;
W3=41984.0/14175.0;
W4=-18160.0/14175.0;
/* INIZIALIZZO LE SOMME CORRENTI A ZERO */
*FLAG=0.0;
*RESULT=0.0;
COR11=0.0;
*ERREST=0.0;
AREA=0.0;
*NOFUN=0;
if(A == B)
return;
/* *** PASSO 2 *** INIZIALIZZAZIONE RELATIVA AL PRIMO INTERVALLO */
LEV=0;
NIM=1;
X0=A;
X[16]=B;
QPREV=0.0;
F0=FUN(X0);
STONE=(B-A)/16.0;
X[8]=(X0+X[16])/2.0;
X[4]=(X0+X[8])/2.0;
X[12]=(X[8]+X[16])/2.0;
X[2]=(X0+X[4])/2.0;
X[6]=(X[4]+X[8])/2.0;
X[10]=(X[8]+X[12])/2.0;
X[14]=(X[12]+X[16])/2.0;
for (J=2; J<=16; J=J+2){
F[J]=FUN(X[J]);
}
*NOFUN=9;
/* *** PASSO 3 *** PARTE CENTRALE RIGUARDANTE I CALCOLI
QUESTI CALCOLI IMPIEGANO I VALORI DI QPREV,X0,X2,X4,...X16,F0,F2,...F16.
SI DETERMINANO I VALORI X1,X3,...X15,F1,F3,...F15,QLEFT,QRIGTH,QNOW,
QDIFF,AREA. */
trenta:
X[1]=(X0+X[2])/2.0;
F[1]=FUN(X[1]);
for (J=3; J<=15; J=J+2){
X[J]=(X[J-1]+X[J+1])/2.0;
F[J]=FUN(X[J]);
}
*NOFUN=*NOFUN+8;
STEP=(X[16]-X0)/16.0;
QLEFT=(W0*(F0+F[8])+W1*(F[1]+F[7])+W2*(F[2]+F[6])+W3*(F[3]+F[5])
+W4*F[4])*STEP;
QRIGHT[LEV+1]=(W0*(F[8]+F[16])+W1*(F[9]+F[15])+W2*(F[10]+F[14])
+W3*(F[11]+F[13])+W4*F[12])*STEP;
QNOW=QLEFT+QRIGHT[LEV+1];
QDIFF=QNOW-QPREV;
AREA=AREA+QDIFF;
/* *** PASSO 4 *** TEST SULLA CONVERGENZA */
ESTERR=fabs(QDIFF)/1023.0;
if (ABSERR > (RELERR*fabs(AREA))*(STEP/STONE))
TOLERR=ABSERR;
else
TOLERR=(RELERR*fabs(AREA))*(STEP/STONE);
if (LEV < LEVMIN)
goto cinquanta;
if (LEV >= LEVMAX)
goto sessantadue;
if (*NOFUN > NOFIN)
goto sessanta;
if (ESTERR <= TOLERR)
goto settanta;
/* *** PASSO 5 *** CASO DELLA NON CONVERGENZA */
cinquanta:
NIM=2*NIM;
LEV=LEV+1;
/* SI MEMORIZZA CIO' CHE PUO' SERVIRE PIU' AVANTI */
for (I=1; I<=8; I++) {
FSAVE[LEV]=F[I+8];
XSAVE[LEV]=X[I+8];
}
/* SI ELABORA CIO' CHE SERVE ORA */
QPREV=QLEFT;
for (I=1; I<=8; I++) {
J=-I;
F[2*J+18]=F[J+9];
X[2*J+18]=X[J+9];
}
goto trenta;
/* *** PASSO 6 *** ISTRUZIONI RELATIVE AL CASO IN CUI CI SONO DELLE
DIFFICOLTA'.
CIOE' QUANDO IL NUMERO DEI VALORI DELLA FUNZIONE USATI PER
DETERMINARE RESULT SUPERA I LIMITI POSTI */
sessanta:
NOFIN=2*NOFIN;
LEVMAX=LEVOUT;
*FLAG=*FLAG+((B-X0)/(B-A));
goto settanta;
/* IL LIVELLO CORRENTE E' LEVMAX */
sessantadue:
*FLAG=*FLAG+1.0;
/* *** PASSO 7 *** CASO DELLA CONVERGENZA
SOMMO I CONTRIBUTI TROVATI ALLE SOMME CORRENTI */
settanta:
*RESULT=*RESULT+QNOW;
*ERREST=*ERREST+ESTERR;
COR11=COR11+QDIFF/1023.0;
/* DETERMINO IL PROSSIMO INTERVALLO */
while (NIM%2 != 0) {
NIM=NIM/2;
LEV=LEV-1;
}
NIM=NIM+1;
if (LEV <= 0)
goto ottanta;
/* ELABORO GLI ELEMENTI RELATIVI AL PROSSIMO INTERVALLO */
QPREV=QRIGHT[LEV];
X0=X[16];
F0=F[16];
for (I=1; I<=8; I++) {
F[2*I]=FSAVE[LEV];
X[2*I]=XSAVE[LEV];
}
goto trenta;
/* *** PASSO 8 *** DETERMINO GLI ULTIMI RISULTATI */
ottanta:
*RESULT=*RESULT+COR11;
if (*ERREST == 0.0)
return;
while (fabs(*RESULT)+(*ERREST) == fabs(*RESULT))
*ERREST=2.0*(*ERREST);
return;
}
#include "quanc8.h"
#include <math.h>
void quanc8(double (*FUN)(double), double A, double B, double ABSERR,
double RELERR, double *RESULT, double *ERREST, int *NOFUN, double *FLAG)
/*
Integrale di FUN(X) da A a B
Routine automatica di tipo adattivo basata
sulle formule di Newton-Cotes ad 8 intervalli.
Input
FUN nome della funzione integranda che viene fornita nella
function FUN(X) di tipo double con argomento double.
A limite inferiore di integrazione.
B limite superiore di integrazione.
RELERR errore relativo.
ABSERR errore assoluto.
Output
RESULT una approssimazione dell'integrale che dovrebbe soddisfare
la precisione data dal meno rigoroso dei due errori
ERREST stima della grandezza dell'errore commesso
NOFUN numero di volte in cui la funzione e' stata calcolata
FLAG indicatore di attendibilita'; se FLAG e' zero si e'
ottenuta la precisione richiesta.
se FLAG e' xxx.yyy allora xxx= il numero di intervalli in cui non e'
stata raggiunta la precisione richiesta, mentre 0.yyy ci da' una
indicazione sul punto nell'intorno del quale non si e' raggiunta la
precisione richiesta.
*/
{
double QRIGHT[32],F[17],X[17],FSAVE[9][31],XSAVE[9][31];
int LEVMIN,LEVMAX,LEVOUT,NOMAX,NOFIN,LEV,NIM,J,I;
double W0,W1,W2,W3,W4,COR11,AREA,X0,F0,STONE,STEP;
double QLEFT,QNOW,QDIFF,QPREV,TOLERR,ESTERR;
/* *** PASSO 1 *** INIZIALIZZAZIONE GENERALE COSTANTI */
LEVMIN=1;
LEVMAX=30;
LEVOUT=6;
NOMAX=5000;
NOFIN=NOMAX-(8*(LEVMAX-LEVOUT+pow(2,(double)(LEVOUT+1))));
/* ISTRUZIONI RIGUARDO AL CASO IN CUI NOFUN RAGGIUNGE IL VALORE NOFIN */
W0=3956.0/14175.0;
W1=23552.0/14175.0;
W2=-3712.0/14175.0;
W3=41984.0/14175.0;
W4=-18160.0/14175.0;
/* INIZIALIZZO LE SOMME CORRENTI A ZERO */
*FLAG=0.0;
*RESULT=0.0;
COR11=0.0;
*ERREST=0.0;
AREA=0.0;
*NOFUN=0;
if(A == B)
return;
/* *** PASSO 2 *** INIZIALIZZAZIONE RELATIVA AL PRIMO INTERVALLO */
LEV=0;
NIM=1;
X0=A;
X[16]=B;
QPREV=0.0;
F0=FUN(X0);
STONE=(B-A)/16.0;
X[8]=(X0+X[16])/2.0;
X[4]=(X0+X[8])/2.0;
X[12]=(X[8]+X[16])/2.0;
X[2]=(X0+X[4])/2.0;
X[6]=(X[4]+X[8])/2.0;
X[10]=(X[8]+X[12])/2.0;
X[14]=(X[12]+X[16])/2.0;
for (J=2; J<=16; J=J+2){
F[J]=FUN(X[J]);
}
*NOFUN=9;
/* *** PASSO 3 *** PARTE CENTRALE RIGUARDANTE I CALCOLI
QUESTI CALCOLI IMPIEGANO I VALORI DI QPREV,X0,X2,X4,...X16,F0,F2,...F16.
SI DETERMINANO I VALORI X1,X3,...X15,F1,F3,...F15,QLEFT,QRIGTH,QNOW,
QDIFF,AREA. */
trenta:
X[1]=(X0+X[2])/2.0;
F[1]=FUN(X[1]);
for (J=3; J<=15; J=J+2){
X[J]=(X[J-1]+X[J+1])/2.0;
F[J]=FUN(X[J]);
}
*NOFUN=*NOFUN+8;
STEP=(X[16]-X0)/16.0;
QLEFT=(W0*(F0+F[8])+W1*(F[1]+F[7])+W2*(F[2]+F[6])+W3*(F[3]+F[5])
+W4*F[4])*STEP;
QRIGHT[LEV+1]=(W0*(F[8]+F[16])+W1*(F[9]+F[15])+W2*(F[10]+F[14])
+W3*(F[11]+F[13])+W4*F[12])*STEP;
QNOW=QLEFT+QRIGHT[LEV+1];
QDIFF=QNOW-QPREV;
AREA=AREA+QDIFF;
/* *** PASSO 4 *** TEST SULLA CONVERGENZA */
ESTERR=fabs(QDIFF)/1023.0;
if (ABSERR > (RELERR*fabs(AREA))*(STEP/STONE))
TOLERR=ABSERR;
else
TOLERR=(RELERR*fabs(AREA))*(STEP/STONE);
if (LEV < LEVMIN)
goto cinquanta;
if (LEV >= LEVMAX)
goto sessantadue;
if (*NOFUN > NOFIN)
goto sessanta;
if (ESTERR <= TOLERR)
goto settanta;
/* *** PASSO 5 *** CASO DELLA NON CONVERGENZA */
cinquanta:
NIM=2*NIM;
LEV=LEV+1;
/* SI MEMORIZZA CIO' CHE PUO' SERVIRE PIU' AVANTI */
for (I=1; I<=8; I++) {
FSAVE[LEV]=F[I+8];
XSAVE[LEV]=X[I+8];
}
/* SI ELABORA CIO' CHE SERVE ORA */
QPREV=QLEFT;
for (I=1; I<=8; I++) {
J=-I;
F[2*J+18]=F[J+9];
X[2*J+18]=X[J+9];
}
goto trenta;
/* *** PASSO 6 *** ISTRUZIONI RELATIVE AL CASO IN CUI CI SONO DELLE
DIFFICOLTA'.
CIOE' QUANDO IL NUMERO DEI VALORI DELLA FUNZIONE USATI PER
DETERMINARE RESULT SUPERA I LIMITI POSTI */
sessanta:
NOFIN=2*NOFIN;
LEVMAX=LEVOUT;
*FLAG=*FLAG+((B-X0)/(B-A));
goto settanta;
/* IL LIVELLO CORRENTE E' LEVMAX */
sessantadue:
*FLAG=*FLAG+1.0;
/* *** PASSO 7 *** CASO DELLA CONVERGENZA
SOMMO I CONTRIBUTI TROVATI ALLE SOMME CORRENTI */
settanta:
*RESULT=*RESULT+QNOW;
*ERREST=*ERREST+ESTERR;
COR11=COR11+QDIFF/1023.0;
/* DETERMINO IL PROSSIMO INTERVALLO */
while (NIM%2 != 0) {
NIM=NIM/2;
LEV=LEV-1;
}
NIM=NIM+1;
if (LEV <= 0)
goto ottanta;
/* ELABORO GLI ELEMENTI RELATIVI AL PROSSIMO INTERVALLO */
QPREV=QRIGHT[LEV];
X0=X[16];
F0=F[16];
for (I=1; I<=8; I++) {
F[2*I]=FSAVE[LEV];
X[2*I]=XSAVE[LEV];
}
goto trenta;
/* *** PASSO 8 *** DETERMINO GLI ULTIMI RISULTATI */
ottanta:
*RESULT=*RESULT+COR11;
if (*ERREST == 0.0)
return;
while (fabs(*RESULT)+(*ERREST) == fabs(*RESULT))
*ERREST=2.0*(*ERREST);
return;
}
[/INDENT]Код программ SPLINE, SEVAL:[INDENT]splines.h
Код:
#ifndef _SPLINES
#define _SPLINES
void spline(int n,double * x,double * y,double * b,double * c,double * d);
double seval(int n, double * u,
double * x,double * y,double * b,double * c,
double * d);
#endif
#define _SPLINES
void spline(int n,double * x,double * y,double * b,double * c,double * d);
double seval(int n, double * u,
double * x,double * y,double * b,double * c,
double * d);
#endif
splines.cpp
Код:
#include "splines.h"
#include <math.h>
void spline(int n,double * x,double * y,double * b,double * c,double * d)
{
int i,ib,nm1;
double t ;
nm1=n-1;
if (n < 2) return;
if (n < 3) goto l20;
d[1]=x[2]-x[1];
c[2]=(y[2]-y[1])/d[1];
for (i=2 ;i<=nm1;i++)
{
d=x[i+1]-x;
b=2*(d[i-1]+d);
c[i+1]=(y[i+1]-y)/d;
c=c[i+1]-c;
}
b[1]=-d[1];
b[n]=-d[n-1];
c[1]=0;
c[n]=0;
if (n == 3) goto l10;
c[1]=c[3]/(x[4]-x[2])-c[2]/(x[3]-x[1]);
c[n]=c[n-1]/(x[n]-x[n-2])-c[n-2]/(x[n-1]-x[n-3]);
c[1]=c[1]*sqrt(d[1])/(x[4]-x[1]);
c[n]=-c[n]*sqrt(d[n-1])/(x[n]-x[n-3]);
l10: for (i=2;i<=n;i++)
{
t=d[i-1]/b[i-1];
b=b-t*d[i-1];
c=c-t*c[i-1];
}
c[n]=c[n]/b[n];
for (ib=1;ib<=nm1;ib++)
{
i=n-ib;
c=(c-d*c[i+1])/b;
}
b[n]=(y[n]-y[nm1])/d[nm1]+d[nm1]*(c[nm1]+2*c[n]);
for (i=1;i<=nm1;i++)
{
b=(y[i+1]-y)/d-d*(c[i+1]+2*c);
d=(c[i+1]-c)/d;
c=3*c;
}
c[n]=3*c[n];
d[n]=d[n-1];
return;
l20: b[1]=(y[2]-y[1])/(x[2]-x[1]);
c[1]=0;
d[1]=0;
b[2]=b[1];
c[2]=0;
d[2]=0;
l30:
return;
}
double seval(int n, double * u,
double * x,double * y,double * b,double * c,
double * d)
{
int i,j,k;
double dx;
i=1;
if (i >= n) i=1;
if (*u < x) goto l10;
if (*u <= x[i+1]) goto l30;
l10: i=1;
j=n+1;
l20: k=(i+j) / 2;
if (*u < x[k]) j=k;
if (*u >= x[k]) i=k;
if (j > (i+1)) goto l20;
l30: dx=*u-x;
return y+dx*(b+dx*(c+dx*d));
}
#include <math.h>
void spline(int n,double * x,double * y,double * b,double * c,double * d)
{
int i,ib,nm1;
double t ;
nm1=n-1;
if (n < 2) return;
if (n < 3) goto l20;
d[1]=x[2]-x[1];
c[2]=(y[2]-y[1])/d[1];
for (i=2 ;i<=nm1;i++)
{
d=x[i+1]-x;
b=2*(d[i-1]+d);
c[i+1]=(y[i+1]-y)/d;
c=c[i+1]-c;
}
b[1]=-d[1];
b[n]=-d[n-1];
c[1]=0;
c[n]=0;
if (n == 3) goto l10;
c[1]=c[3]/(x[4]-x[2])-c[2]/(x[3]-x[1]);
c[n]=c[n-1]/(x[n]-x[n-2])-c[n-2]/(x[n-1]-x[n-3]);
c[1]=c[1]*sqrt(d[1])/(x[4]-x[1]);
c[n]=-c[n]*sqrt(d[n-1])/(x[n]-x[n-3]);
l10: for (i=2;i<=n;i++)
{
t=d[i-1]/b[i-1];
b=b-t*d[i-1];
c=c-t*c[i-1];
}
c[n]=c[n]/b[n];
for (ib=1;ib<=nm1;ib++)
{
i=n-ib;
c=(c-d*c[i+1])/b;
}
b[n]=(y[n]-y[nm1])/d[nm1]+d[nm1]*(c[nm1]+2*c[n]);
for (i=1;i<=nm1;i++)
{
b=(y[i+1]-y)/d-d*(c[i+1]+2*c);
d=(c[i+1]-c)/d;
c=3*c;
}
c[n]=3*c[n];
d[n]=d[n-1];
return;
l20: b[1]=(y[2]-y[1])/(x[2]-x[1]);
c[1]=0;
d[1]=0;
b[2]=b[1];
c[2]=0;
d[2]=0;
l30:
return;
}
double seval(int n, double * u,
double * x,double * y,double * b,double * c,
double * d)
{
int i,j,k;
double dx;
i=1;
if (i >= n) i=1;
if (*u < x) goto l10;
if (*u <= x[i+1]) goto l30;
l10: i=1;
j=n+1;
l20: k=(i+j) / 2;
if (*u < x[k]) j=k;
if (*u >= x[k]) i=k;
if (j > (i+1)) goto l20;
l30: dx=*u-x;
return y+dx*(b+dx*(c+dx*d));
}
Среда разработки: Visual Studio 2008
Операционная система: Windows 7, а также проверенно на Windows XP
В опреациях деления, проверь чтобы делитель не являля 0.0.
Думаю поможет.
Думаю поможет.
Деление на ноль отметается...
Эти реализации кванка/севала/сплайна проверены в течении нескольких лет. Этот код выдавал преподаватель. Файлы несколько раз перезаливал (вдруг я накасячил когда код листал) не помогает.
Эти реализации кванка/севала/сплайна проверены в течении нескольких лет. Этот код выдавал преподаватель. Файлы несколько раз перезаливал (вдруг я накасячил когда код листал) не помогает.
Разбирайся в арифметике... где то получается минус бесконечность.
DragonHT, а программная функция
представляет ту самую математическую функцию, решение для которой вы хотите найти? Если да, то как-то странным выглядит использование в ней, для вычисления результата, глобального массива x, и глобальной переменной-индекса i. Как мне думается, в подобных случаях возвращаемый программной функцией результат должен целиком и полностью зависеть только от передаваемых в функцию аргументов (в данном случае - от double t). Так что, может быть, надо умножение на x удалить, оставить только exp(t)*sin(t).
И еще учтите, что здесь у вас используется экспонента от x (функция exp из библиотеки math), которая очень быстро возрастает при увеличении аргумента. Возможно из-за этого у вас и возникает то о чем говорят предыдущие авторы (то что где-то получается бесконечность).
Посмотрел также на код функции seval, на которую вы грешите. Прежде всего бросается в глаза очевидная нелепость:
Вторая операция не имеет здесь абсолютно никакого смысла. Почему - предлагаю вам подумать самостоятельно.
Есть также замечания по структурированию программного кода в следующем фрагменте:
Во-первых, использование операторов goto считается дурным тоном, здесь вполне можно обойтись без goto.
Во-вторых, в блоке обозначенном меткой 120, есть два условия if (*u < x[k]) и (*u >= x[k]). Так вот второе можно не писать, а написать else после первого, т.к. если выполняется первое, то не выполняется второе, и наоборот.
Вот полный эквивалент этого фрагмента кода, с устранением этих замечаний:
Думаю, не будете возражать, что если так переписать, то логика работы программы становится гораздо более понятней чем в том варианте, что у вас? :)
Код:
double FUN ( double t)
{
return exp(x*t)*sin(t);
}
{
return exp(x*t)*sin(t);
}
представляет ту самую математическую функцию, решение для которой вы хотите найти? Если да, то как-то странным выглядит использование в ней, для вычисления результата, глобального массива x, и глобальной переменной-индекса i. Как мне думается, в подобных случаях возвращаемый программной функцией результат должен целиком и полностью зависеть только от передаваемых в функцию аргументов (в данном случае - от double t). Так что, может быть, надо умножение на x удалить, оставить только exp(t)*sin(t).
И еще учтите, что здесь у вас используется экспонента от x (функция exp из библиотеки math), которая очень быстро возрастает при увеличении аргумента. Возможно из-за этого у вас и возникает то о чем говорят предыдущие авторы (то что где-то получается бесконечность).
Посмотрел также на код функции seval, на которую вы грешите. Прежде всего бросается в глаза очевидная нелепость:
Код:
i=1;
if (i >= n) i=1;
if (i >= n) i=1;
Вторая операция не имеет здесь абсолютно никакого смысла. Почему - предлагаю вам подумать самостоятельно.
Есть также замечания по структурированию программного кода в следующем фрагменте:
Код:
if (*u < x) goto l10;
if (*u <= x[i+1]) goto l30;
l10: i=1;
j=n+1;
l20: k=(i+j) / 2;
if (*u < x[k]) j=k;
if (*u >= x[k]) i=k;
if (j > (i+1)) goto l20;
l30: dx=*u-x;
return y+dx*(b+dx*(c+dx*d));
if (*u <= x[i+1]) goto l30;
l10: i=1;
j=n+1;
l20: k=(i+j) / 2;
if (*u < x[k]) j=k;
if (*u >= x[k]) i=k;
if (j > (i+1)) goto l20;
l30: dx=*u-x;
return y+dx*(b+dx*(c+dx*d));
Во-первых, использование операторов goto считается дурным тоном, здесь вполне можно обойтись без goto.
Во-вторых, в блоке обозначенном меткой 120, есть два условия if (*u < x[k]) и (*u >= x[k]). Так вот второе можно не писать, а написать else после первого, т.к. если выполняется первое, то не выполняется второе, и наоборот.
Вот полный эквивалент этого фрагмента кода, с устранением этих замечаний:
Код:
if(*u < x || *u > x[i+1])
{
i = 1;
j = n + 1;
do
{
k = (i + j) / 2;
if (*u < x[k])
j = k;
else
i = k;
}
while(j > (i + 1));
}
dx = *u - x;
return y + dx * (b + dx * (c + dx * d));
{
i = 1;
j = n + 1;
do
{
k = (i + j) / 2;
if (*u < x[k])
j = k;
else
i = k;
}
while(j > (i + 1));
}
dx = *u - x;
return y + dx * (b + dx * (c + dx * d));
Думаю, не будете возражать, что если так переписать, то логика работы программы становится гораздо более понятней чем в том варианте, что у вас? :)