Метод конечных разностей

[manyak666]

Здравствуйте! Помогите пожалуйста! Надо написать прогу на Си, которая решает краевую задачу для обыкновенного дифференциального уравнения 2 порядка методом конечных разностей. Никак не могу разобраться уже неделю( Подскажите хотя бы, с чего начать и как вообще её сделать? Или, если у кого-то уже есть такая, поделитесь пожалуйста!

[manyak666]

Ну что, никто не может подсказать чтоли?(

[Abstraction]

Можете заглянуть в эту тему.

[Bugrimov]

Abstraction подкинул неплохую тему. А вообще нужно начинать с решения задачи на листе бумаги. Реализация на языке программирования - это уже дело техники.

[Abstraction]

А, да. В фейнмановких лекциях где-то был хороший пример.
Вот, нашлось.

[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();
}