Метод конечных разностей для краевой задачи

[Manyak666]

Здравствуйте! Кому-нибудь приходилось решать краевую задачу для дифференциального уравнения методом конечных разностей (он же метод сеток)? Если есть готовая программа, поделитесь пожалуйста, буду очень признателен) Помогите, уже крышак едет(

[Manyak666]

В общем нашлась готовая программа на делфи в виде приложения с интерфейсом. Вместе с другом-программистом перепилили её на С++ в борланде, после двух недель страдания она заработала. Теперь нужно её переделать в обычное консольное приложение, и чтобы работало в Visual C++, ибо разбираться с борландом - это для меня ещё одна трудно преодолимая проблема (борланд стоял у коллеги). В итоге у нас чё-то получилось, но оно коряво работает и выдаёт полную ерунду. Помогите, пожалуйста, разобраться и исправить ошибки.
Вот наше "творчество":

Код: Выделить всё

#pragma hdrstop
#include <stdio.h>
#include <conio.h>

const int N = 10;
const int M = 30;
double eps = 0.001;
double F0[N][M], F[N][M], Xl[N];

void exchange(double **mas, int razmer,int from, int to){
for(int i=0; i < razmer; i++){
mas[from][i] = mas[from][i] + mas[to][i];
mas[to][i] = mas[from][i] - mas[to][i];
mas[from][i] = mas[from][i] - mas[to][i];
}
}

// определяем ближайшую строку с не 0 диоганалью вниз
int get_n_no_empty(double **a, int razmer, int begin){
for(int i=begin+1; i < razmer; i++){
if(a[i][i] != 0){
return i;
}
}
}

double *Gauss(double **a, double *b)
{
int razmer; //размер матрицы
double *x;

/* Здесь переменной razmer необходимо присвоить размер матрицы */
// создание динамических переменных
/* Здесь нам надо присвоить значения матрице а, и заполнить свободные члены b */

// прямой проход
// идем слево на право по коэфицентам х
for(int k=0; k < razmer-1; k++){
// вычитаем по строкам
for(int m=k+1; m < razmer; m++){
// если на диагонали элемент = 0, то поменяем местами строки
if(a[m][m] == 0) exchange(a, razmer, m, get_n_no_empty(a, razmer, m));

double koeficent=a[m][k]/a[0][k];
//вычисление новых коэфицентов уравнения
b[m] = b[m] - b[0] * koeficent;
for(int z=0; z < razmer; z++){
a[m][z] = a[m][k] - a[0][k] * koeficent;
}

}
}
// ищем решения
for(int m=razmer-1; m >= 0; m--){
double sum=0;
// идем по строке спаво налево, считая сумму корень*коэфицент, до текущего корня
for(int i=razmer-1; i > m; i--){
sum += x[i] * a[m][i];
}
x[m] = (b[m] - sum)/a[m][m];
}

// вывод решений
// cout « "###Решения:" « endl;
return x;
}

double Fn(double t)
{
double x = t*t-1;
return x;
}

void progon(double E[N][M], double D[N])
{
int i;
// double znam;
double a[N],b[N],c[N],f[N],alpha[N],beta[N];
// int nn;
double X[N];
// double G[N][M];

for (i=0;i<=N;i++)
{
a[i] = E[i+1][i];
c[i] = E[i][i];
b[i] = E[i-1][i];
}
for (i=0;i<=N;i++)
{
f[i] = D[i];
}
c[0] = 1;
c[N] = 1;
alpha[0] = 1;
beta[0] = 0;
for (i=1;i<=N-1;i++)
{
alpha[i+1] = b[i]/(c[i]-a[i]*alpha[i]);
beta[i+1] = (a[i]*beta[i] + f[i])/(c[i]-a[i]*alpha[i]);
}
for (i=N-1;i>=0;i--)
{
X[i] = alpha[i+1]*X[i+1] + beta[i+1];
}
for (i=0;i<=N;i++)
{
Xl[i] = X[i];
}
}

void main(){

int i,j,k;
double h,tau;
double a,b,Xi,Tj,a1;
double E[N][M];
double D[N],*Ms;

a = 0; b = 1;
h = (b-a)/N;
tau = h/2;
for (i=1;i<=M;i++)
{
F[0][i] = 0.01;
F[N][i] = 0;
}
for (j=0;j<=N;j++)
{
F[j][0] = 0.01*(1-j*h);
}

for (j=1;j<=M;j++)
{
E[0][0] = 0;
E[0][1] = 1;
E[N][N] = 1;
E[N-1][N] = 0;
D[0] = 0.01;
D[N] = 0;
for (i=1;i<=N-1;i++)
{
Xi = i*h;
Tj = tau*j;
a1 = 3.3-1.5*Xi;
E[i-1][i] = -(a1*a1)/(2*h*h);
E[i][i] = 1/tau + (a1*a1)/(h*h);
E[i+1][i] = -(a1*a1)/(2*h*h);
D[i] = Fn(Tj) + (a1*a1)/(2*h*h)*F[i+1][j] + (1/tau - (a1*a1)/(h*h))*F[i][j] + (a1*a1)/(2*h*h)*F[i-1][j];
}
// Ms = Gauss(E,D);
progon(E,D);
for (k = 1;k<=N-1;k++)
{
F[k][j] = Xl[k];
}

for (i=0;i<=N;i++)
{
for (j=0;j<=M;j++)
{
printf("\n%f", F[i][j]);

}
}

}
getch();
}

[Romeo]

Дело в том, что я теории всей не помню, так что формулы проверять не возьмусь.

А могло так получится, что при переносе из Борланда в консольное приложение просто-напросто где-то ошибку в формуле допустили? Если есть один работающий код, и есть второй переписанный, но не работающий, это ведь совсем не проблема найти ошибку. Просто внимательно проверить правильность переноса и всё.

[Manyak666]

Дело в том, что у меня нет борланда и, к тому же, я не совсем всё понял, как мы это переносили, а коллега-программист уехал куда-то отдыхать, и компьютера у него там нет.. В формулах ошибок не должно быть, мы основные функции не трогали, они вроде как везде одинаково выглядят

[Manyak666]

Сделали другую программу, только она какую-то странную погрешность выдаёт, посмотрите, пожалуйста, что не так..

Код: Выделить всё

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
using namespace std;

void main()
{
	int n,i,k;
	float h,c1,c2,c,d1,d2,d,a,b,H,eps;
	float p,q, fx[9999], betta[9999],gamma[9999],fi[9999],alfa[9999];
	float y[9999],u[9999],v[9999],x[9999],Y[9999];
	c1=1.0;
	c2=-1.0;
	d1=1.0;
	d2=2.71828183;
	a=0; b=1;
	n=10;
	c=0; d=1;
	x[0]=a;
	h=(b-a)/n;
	betta[0]=c1*h-c2;
	gamma[0]=c2;
	fi[0]=h*c;
	alfa[0]=0;
	for(i=1;i<n;i++)
	{
		x[i]=x[i-1]+h;
		p=1;
		q=-1;
		fx[i] = -exp(-x[i])*cos(x[i]);
		fi[i]=fx[i]*h*h;
		alfa[i]=1-(1/2)*p;
		betta[i]=q*h*h-2;
		gamma[i]=1+(1/2)*p*h;
	}
	alfa[n]=-d2;
	betta[n]=h*d1+d2;
	fi[n]=h*d;
	v[0]=gamma[0]/betta[0];
	u[0]=fi[0]/betta[0];
	for(i=1;i<n;i++)
	{
		v[i]=-gamma[i]/(betta[i]+alfa[i]*v[i-1]);
		u[i]=(fi[i]-alfa[i]*u[i-1])/(betta[i]+alfa[i]*v[i-1]);
	}
	y[n]=u[n];
	for(i=n-1;i>0;i--) 
		y[i]=u[i]+v[i]*y[i+1];
	cout<<"   i"<<"  x  "<<"     y"<<endl;
	for(i=0;i<n;i++)
	{
		cout.width(4);
		cout<<i;
		cout.width(5);
		cout<<x[i];
		cout.width(14);
		cout<<y[i]<<endl;
	}

	H=(b-a)/(2*n);
	betta[0]=c1*H-c2;
	gamma[0]=c2;
	fi[0]=H*c;
	alfa[0]=0;
	k=2*n;
	for(i=1;i<k;i++)
	{
		x[i]=x[i-1]+H;
		p=1;
		q=-1;
		fx[i] = -exp(-x[i])*cos(x[i]);
		fi[i]=fx[i]*H*H;
		alfa[i]=1-(1/2)*p;
		betta[i]=q*H*H-2;
		gamma[i]=1+(1/2)*p*H;
	}
	alfa[k]=-d2;
	betta[k]=H*d1+d2;
	fi[k]=H*d;
	v[0]=gamma[0]/betta[0];
	u[0]=fi[0]/betta[0];
	for(i=1;i<k;i++)
	{
		v[i]=-gamma[i]/(betta[i]+alfa[i]*v[i-1]);
		u[i]=(fi[i]-alfa[i]*u[i-1])/(betta[i]+alfa[i]*v[i-1]);
	}
	Y[k]=u[k];
	for(i=2*n-1;i>0;i--) 
		Y[i]=u[i]+v[i]*Y[i+1];
	cout<<"   i"<<"  x  "<<"     Y"<<endl;
	for(i=0;i<k;i++)
	{
		cout.width(4);
		cout<<i;
		cout.width(5);
		cout<<x[i];
		cout.width(14);
		cout<<Y[i]<<endl;
	}
	cout<<"eps:"<<endl;
	for(i=0;i<n;i++)
	{
		eps=(Y[2*i]-y[i]);
		cout<<eps<<endl;
	}
	system("PAUSE");
}

[somewhere]

Может тебе проще еще и исходник на дельфе приложить? А то бывали случаи, при переносе, кодеры не правильно трактовали синтаксис языка и ничего потом не работало. Особенно если в программе были указатели.

[Manyak666]

Лучше помогите со второй, которую мы делали сами полностью, пожалуйста! С делфи я совсем не знаком, а в С++ уже чё-то понимаю..