Re: Кто напишет прогу по ОММ за вознаграждение?

mamycik

Очень нужно сделать две задачи по ОММ для экзамена (3 курс физфак)! Экзамен в среду. Вознаграждение гарантируется! тел. 89265825410, звонить в любое время!

SvinkaVJeansah

Обрисуй, что за прога в двух словах.

garikus

Какой номер у второго задания ?

mamycik

Задание 2
№106
Используя метод конечных разностей (неявную схему решить краевую задачу:
du/dt=d2u/dx2+cosx exp(-t x: от 0 до п, t больше 0
при х=0 u=0
при х=п u=0
при t=0 u=cosx
прога должна реализовывать метод прогонки.
а еще надо решить численно задачу 1(г Гл. 3 "Задач по мат. физике" Боголюбова, Кравцова
Помогите плиз, очень надо!!

SvinkaVJeansah

Могу написать то заданее, что выложено здесь. Что там в книге Боголюбова - не знаю.

SvinkaVJeansah

89265752211 )

mamycik

Супер! напиши тогда свой телефон или комнату, я зайду-договоримся

garikus

если что - можешь ко мне зайти - прога с графиками готова
Можно язык выбрать: C; Turbo Pascal / Delphi; Oberon-2 / Component Pascal

Bird_V

Почему-то мне помнится, что 1е задание (из книги Боголюбова) надо было в МатЛабе () делать в пакете каком-то () для решения Дифференциальных Уравнений в Частных Производных...

al70

Для начала нужно согласовать начальные условия с граничными. Ты или в начальном условии cos(x) на sin(x) поменяй или граничные в 1 и -1 выстави.

mamycik

а ты сейчас не спишь еще?

SvinkaVJeansah



// ***************************************************

// 

//	Ut = Uxx + cos(x)*exp(-t); 0<=x<=pi; t>0; 

//	U(t,0) = cos(x);

//	U(0,x) = 0;

//	Ut(0,x) = 0;

//

// ***************************************************

#include <math.h>

#include <iostream.h>



#define NT	120

#define	NX	40

#define pi	3.1459



// for single-vector calc

double a[NX], b[NX], u0[NX], u1[NX];

// array

double u[NT][NX];

double maxx = pi;

double maxt = 10;





void calc(double, double, double);

void print(double, double);



void main(void)

{

	int i, j;

	double dt = maxt/(NT-1);

	double dx = maxx/(NX-1);

	double t=0;

		

	// start calc

	for(i = 0; i < NX; i++)

		u[0][i] = cos(dx * i)*exp(-0);

		

	// steb by step: t

	for(i = 1; i < NT; i++)

	{

		t+=dt;



		// init vectors

		for(j = 0; j < NX; j++)

			u0[j] = u[i-1][j];



		u1[0] = 0;

		u1[NX - 1] = 0;

		

		// calc

		calc(dt, dx, t);

		

		// come back

		for(j = 0; j < NX; j++)

			u[i][j] = u1[j];

	}



		print(dt, dx);

}



void calc(double dt, double dx, double t)

{

	int i;



	a[1] = dt/(dx*dx + 2*dt); 

	b[1] = cos(dx)*exp(-t) + u1[0]/(dx*dx) + u0[1]/dt;

	b[1] = b[1]*(dt*dx*dx)/(dx*dx + 2*dt);



	for(i = 2; i < NX - 1; i++)

	{

		a[i] = dt/(dx*dx + 2*dt - a[i-1]*dt);

		b[i] = cos(dx*i)*exp(-t) + u0[i]/dt + b[i-1]/(dx*dx);

		b[i] = b[i]*dt*dx*dx/(dx*dx + 2*dt - a[i-1]*dt);

	}

	// back

	for(i = NX-2; i > 0; i--)

		u1[i] = a[i]*u1[i+1] + b[i];



 

}



void print(double dt, double dx)

{

	cout << "(t,x,u)" << endl;

	for(int i=0; i < NT; i++)

		for(int j=0; j<NX; j++)

			cout << i*dt << ", " << j*dx << ", " << u[i][j] << endl;



}

mamycik

Супер, объяснишь тогда завтра, ок?

SvinkaVJeansah

Ага. А почему решил, что супер?!

garikus

По крайней мере график получается таким, каким должен быть (только всё-таки надо заменить cos(x) на sin(x) в условии u(x)|t=0)

SvinkaVJeansah

Типа для согласования граничных\начальных? Можно, хотя один хрен в модели с тремя точками угловые значения не используются

garikus



MODULE Omm;



	(*

		2005 Alexander Shiryaev, 309 aka Aix-D

		Programming language: Component Pascal

		Last modified: 31.05.2005

	*)



	IMPORT Math, TextMappers, TextModels, TextViews, Views;



	CONST

		N = 100;

		M = 100;

		time = 10;



	VAR

		y: ARRAY N + 1, M + 1 OF REAL;

		h, tau: REAL;



	PROCEDURE Calc;

		VAR n, m: INTEGER;

			An, Bn, Cn, Fn: REAL;

			alpha, beta: ARRAY N + 1 OF REAL;



	BEGIN

	

		h := Math.Pi / N;

		tau := time / M;



		FOR n := 0 TO N DO

			y[n][0] := Math.Sin(n * h)

		END;



		FOR m := 0 TO M DO

			y[0][m] := 0; y[N][m] := 0

		END;



		(* Метод прогонки *)



		An := 1 / (h * h);

		Bn := 1 / (h * h);

		Cn := 1 / tau + 2 / (h * h);



		FOR m := 1 TO M DO

		

			(* Прямой ход *)



			alpha[1] := 0; beta[1] := 0;



			FOR n := 1 TO N - 1 DO

				alpha[n + 1] := Bn / (Cn - alpha[n] * An);

				Fn := 1 / tau * (y[n][m - 1]) + Math.Cos(h * n) * Math.Exp(-tau * (m - 1;

				beta[n + 1] := (An * beta[n] + Fn) / (Cn - alpha[n] * An)

			END;



			y[N - 1][m] := 0;



			(* Обратный ход *)



			n := N - 2;

			WHILE n > 0 DO

				y[n][m] := alpha[n + 1] * y[n + 1][m] + beta[n + 1];

				DEC(n)

			END

			

		END

		

	END Calc;



	PROCEDURE Output;

		VAR f: TextMappers.Formatter;

			m: TextModels.Model;

			v: TextViews.View;

			i, j: INTEGER;

	BEGIN

		m := TextModels.dir.New;

		f.ConnectTo(m);



		FOR j := 0 TO M DO

			FOR i := 0 TO N DO

				f.WriteReal(i * h); f.WriteTab;

				f.WriteReal(j * tau); f.WriteTab;

				f.WriteReal(y[i][j]); f.WriteLn

			END

		END;



		v := TextViews.dir.New(m);

		Views.OpenView(v)

	END Output;



	PROCEDURE Do*;

	BEGIN

		Calc; Output

	END Do;



END Omm.

SvinkaVJeansah

перекрестная прогонка )


//******************************************************

//

//	Utt = Uxx + Uyy;

//	0 <= x,y,t <= 1;

//	u(0,y,t) = t*y^2;

//	u(1,y,t) = t(1 + t*y^2);

//	u(x,0,t) = t*x^2;

//	u(x,1,t) = t(1 + t*x^2);

//

//******************************************************



#include <iostream.h>





#define NZ	25

#define NT	25



// for single-vector cflc

double a[NZ], b[NZ], u0[NZ], u1[NZ];

// layer

double u[NZ][NZ];



void calc(double delta_coord, double delta_time);

void print_layer(double time, double delta_x, double delta_y);



void main(void)

{

	double maxx = 1, maxy = 1, maxt = 1;

	double dx = maxx/(NZ - 1 dy = maxy/(NZ - 1 dt = maxt/(NT - 1);



	int i,j;



	// poluting (x,y,0) -> u

	for(i = 0; i < NZ; i++)

			for(j = 0; j < NZ; j++)

				u[i][j] = 0;

				

	int nt = 0;

	bool isXlayer = true;

	int nx, ny;



	while(nt < NT)

	{

		// printing layer out

		print_layer(nt*dt, dx, dy);

		// choosing diredtion

		if(isXlayer)

		{

			isXlayer = false;		

			for(nx = 0; nx < NZ; nx++)

			{

				// layer -> temp

				for(i = 0; i < NZ; i++)

					u0[i] = u[nx][i];

				// border values

				u1[0] = nt * dt * (nx * dx) * (nx * dx);

				u1[NZ-1] = nt * dt * (1 + (nx * dx) * (nx * dx;

				// calc temp

					calc(dy, dt);

				// temp -> layer

				for(i = 0; i < NZ; i++)

					u[nx][i] = u1[i];

			}

		}

		else

		{

			isXlayer = true;

			for(ny = 0; ny < NZ; ny++)

			{

				// layer -> temp

				for(i = 0; i < NZ; i++)

					u0[i] = u[i][ny];

				// border values

				u1[0] = nt * dt * (ny * dy) * (ny * dy);

				u1[NZ-1] = nt * dt * (1 + (ny * dy) * (ny * dy;

				// calc temp

					calc(dx, dt);

				// temp -> layer

				for(i = 0; i < NZ; i++)

					u[i][ny] = u1[i];

			}

		}

	nt++;

	}

}



// calculating vector

void calc(double d, double dt)

{

	// start coef.

	a[1] = dt / (2 * (d*d + dt;

	b[1] = (u1[0]*dt/2 + u0[1]*d*d)/(d*d + dt);



	int i;

	// calc coef.

	for(i = 2; i < NZ-1; i++)

	{

		a[i] = dt / (2*(d*d + dt) - a[i-1]*dt);

		b[i] = (b[i-1]*dt + 2*u0[i]*d*d) / (2*(d*d + dt) - a[i-1]*dt);

	}

	// calc values

	for(i = NZ-2; i > 0; i--)

		u1[i] = a[i]*u1[i+1] + b[i];

}



//printing u[][] layer

void print_layer(double time, double dx, double dy)

{

	for(int i = 0; i < NZ; i++)

		for(int j = 0; j < NZ; j++)

			cout << time << ", " << dx*i << ", " << dy*j << ", " << u[i][j] << endl;

}

mamycik

Осталась задача 1(г) гл.3 боголюбов-Кравцов (численным методом). Очень срочно надо сегодня! Кто-нить сделает за 150-200р.?
Оставить комментарий
Имя или ник:
Комментарий: