Приближённая минимизация интегрального функционала по методу Ритца

Размещено на http://www.allbest.ru/

- 24 -

Кафедра «Прикладная математика»

Курсовая работа

по дисциплине: «Вариационное исчисление»

на тему: «Приближённая минимизация интегрального функционала по методу Ритца»

Москва 2005 г.

Содержание

Постановка задачи

I. Точное решение

1. Основное условие экстремума

2. Решение поставленной задачи

II. Минимизация заданного функционала по методу Ритца

1. Описание метода Ритца

2. Пример применения метода Ритца

3. Решение поставленной задачи - непосредственное отыскание первой итерации

4. Отыскание последовательности итераций по методу Ритца (n = 1,2,3,4,5) с помощью Mathcad 2001

5. Невязки. Графики приближений

6. Результаты работы

III. Приближённое решение краевой задачи для уравнения Эйлера методом конечных разностей

1. Описание метода конечных разностей

2. Составление конечноразностного уравнения

3. Описание метода прогонки

4. Листинг программы на языке С++

5. Расчёт приближений

6. Оценка невязок

Выводы

Список литературы

Постановка задачи

Найти точное решение задачи о минимуме заданного функционала J[y] - экстремаль y0(x). , y(0) = 0, y(1) = sh1

Найти n приближённых решений (итераций) yn(x) задачи о минимуме J[y] по методу Ритца при указанном в задании выборе системы координатных функций цk(x), k = 1, 2, …, n.

Исходная точность 0.00001; далее - по ходу задачи.

Сравнить yn(x) и y0(x) в 10 точках

xi = a + 0.1(b-a)i; i = 0, 1, 2, …, 10,

указав (по 10 точкам) невязки

дn = max| yn(xi) - y0(xi) |, n = 1,2,3,4,5;

Дn = max| yn+1(xi) - yn(xi) |, n = 1,2,3,4.

Убедиться в том, что дn > 0 ( и Дn > 0) с ростом n.

Построить графики yn(x) и y0(x).

Решить приближённо краевую задачу для уравнения Эйлера (рассматривавшегося в п. 1 ) методом конечных разностей ( число шагов последовательно равно n = 10, 20, 100); сравнить приближённое решение с точным, указав невязки - максимум модуля разности приближённого и точного решений.

I. Точное решение

1. Основное условие экстремума

Вариационное исчисление рассматривает решения задач, связанных с минимизацией следующего функционала (а так же с обобщениями этой задачи):

.

Здесь - некоторая непрерывная функция трех переменных; допустимыми являются непрерывно-дифференцируемые функции y(x), удовлетворяющие краевым условиям и .

Имеется основное необходимое условие экстремума основной задачи классического вариационного исчисления, которое в известной мере подобно уравнению в основной задаче математического анализа. Здесь основное необходимое условие экстремума заключается в следующем: функция , дающая минимум (максимум) заданному функционалу, должна удовлетворять следующему дифференциальному уравнению второго порядка, так называемому уравнению Эйлера:

.

Функция так же должна удовлетворять краевым условиям и .

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

2. Решение поставленной задачи

Для заданного функционала составим и решим уравнение Эйлера.

; ; ; ;

;

;

Получилось однородное дифференциальное уравнение второго порядка.

Составим характеристическое уравнение: л2 - 1 = 0.

Его корни

Следовательно решением будет

y = C1 ex+ C2e-x=D1shx+D2chx;

Используем краевые условия

y(0) = 0; D1sh(0)+D2ch(0) > D2 = 0;

y(1) = sh(1); D1sh(1)=sh(1) > D1 = 1;

В итоге точным решением задачи (т. е. функцией, дающей экстремум J[y]) является

y(x) = sh(x)

II. Минимизация заданного функционала по методу Ритца

1. Описание метода Ритца

Основная идея метода состоит в том, что при разыскании экстремума функционала J[y(x)] рассматривается не всё пространство допустимых функций, а лишь всевозможные линейные комбинации допустимых функций вида

(1)

где бi - постоянные, а система {цi(x)}, называемая системой координатных функций, такова, что функции цi(x) линейно независимы и образуют в рассматриваемом пространстве полную систему функций.

Требование, чтобы yn(x) были допустимыми функциями, вообще говоря, накладывает на координатные цi(x) некоторые дополнительные условия типа условий гладкости или удовлетворения граничным условиям. На таких линейных комбинациях функционал J[y(x)] обращается в функцию аргументов б1, б2, б3, …, бn:

J[yn(x)] = Ф(б1, б2, б3, …, бn) (2)

Находим те значения б1, б2, б3, …, бn, которые доставляют функции Ф(б1, б2, б3, …, бn) экстремум; для этого решаем систему, вообще говоря нелинейных относительно б1, б2, б3, …, бn уравнений

(i = 1, 2, …, n) (3)

и найденные значения бi подставляем в (1). Полученная таким образом последовательность {yn(x)} является минимизирующей последовательностью, т. е. такой, для которой последовательность значений функционала { J[yn(x)]} сходится к минимуму или к нижней грани значений функционала J[y(x)]. Однако из того, что

,

ещё не следует, что . Минимизирующая последовательность может и не стремиться к функции, реализующей экстремум в классе допустимых функций.

Можно указать условия, обеспечивающие существование абсолютного минимума функционала и его достижение на функциях {yn(x)}.

В случае, когда ищется экстремум функционала

,

y(x1)=y1, y(x2)=y2,

эти условия таковы:

Функция F(x,y,z) непрерывна по совокупности своих аргументов при любом z и при (x,y) Є D, где D - замкнутая область плоскости XOY, в которой лежат линии yn(x)

Существуют константы б > 0, р > 1, в, для которых

,

каково бы ни было z и для любой точки (x,y) Є D.

Функция F(x,y,z) имеет непрерывную частную производную Fz(x,y,z), причём эта производная для любой точки (x,y) Є D есть неубывающая функция от z .

Сформулированные выше условия выполняются, в частности, для функционалов

вида

,

y(x1)=a, y(x2)=b,

где p(x), q(x), r(x) - заданные непрерывные на [x1,x2] функции, причём p(x) имеет непрерывную производную и p(x) > 0, q(x) ? 0.

Если таким методом определяется абсолютный экстремум функционала, то приближённое значение минимума функционала получается с избытком, а максимума - с недостатком. От удачного выбора системы координатных функций {цi(x)} в значительной степени зависит успех применения этого метода.

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

В случае, когда приходится находить приближённо экстремум функционалов J[z(x1,x2, …, xn)], зависящих от функций нескольких независимых переменных, выбирается координатная система функций ц1(x1,x2, …, xn), ц2(x1,x2, …, xn), …, цm(x1,x2, …, xn) и приближённое решение вариационной задачи ищется в виде

,

где коэффициенты бk - некоторые постоянные числа. Для определения их аналогично предыдущему составляем систему уравнений (k = 1, 2, …, m), где Ф(б1, б2, б3, …, бn) - результат подстановки zm в функционал J[z].

2. Пример применения метода Ритца

Требуется найти приближённое решение задачи о минимуме функционала

, (1)

y(0) = 0, y(1) = 0 (2)

и сравнить с точным решением.

Пусть (k = 1, 2, …) - система координатных функций.

Функции цk(x), очевидно, удовлетворяют краевым условиям цk(0) = 0, цk(1) = 0, являются линейно независимыми и представляют в пространстве C1[0,1] полную систему.

При k = 1 получаем y1(x) = б1(x-x2). Подставляя это выражение для y1(x) в функционал (1), получим



Коэффициент б1 находим из уравнения

,

откуда . Следовательно,

.

Рассмотрим точное решение. Уравнение Эйлера для данного функционала:

.

Решая это неоднородное линейное уравнение, находим

.

Используя граничные условия y(0) = 0, y(1) = 0, получим окончательно:

.

В результате сравнения точного решения с приближённым получаем:

x	Точное решение	Приближённое решение
0,00	0	0
0,25	-0,044	-0,052
0,50	-0,070	-0,069
0,75	-0,060	-0,052
1,00	0	0

3.Решение поставленной задачи - непосредственное отыскание первой итерации

, y(0) = 0, y(1) =sh1

ц0(x) = xsh1

ц0(x) = x(1-x)k

y1 = xsh1 + б1(x-x2)

y1(x) = 1,175201x - 0,267091x(1 - x)

X	y0	y1	| y1- y0|
0	0	0	0
0,1	0,100167	0,093482	0,006685
0,2	0,201336	0,192306	0,00903
0,3	0,30452	0,296471	0,008049
0,4	0,410752	0,405979	0,004774
0,5	0,521095	0,520828	0,000268
0,6	0,636654	0,641019	0,004365
0,7	0,758584	0,766552	0,007968
0,8	0,888106	0,897426	0,00932
0,9	1,026517	1,033643	0,007126
1	1,175201	1,175201	0

max|y1- y0| = 0,000932

4. Отыскание последовательности итераций по методу Ритца (n = 1,2,3,4,5) с помощью Mathcad 2001



b1=-0.362747 b2=0.191312



c1=-0.367179 c2=0.213327 c3=-0.22015

d1=-0.3678 d2=0.219442 d3=-0.036272 d4=9.504139x10-3

e1=0.040634 e2=-0.019294 e3=-1.012313x10-3 e4=2.408927x10-4 e5=50356116x10-6



5.Невязки. Графики приближений

Рассмотрим невязки:

(Как говорилось ранее, невязки выражаются следующим образом:

дn = max| yn(xi) - y0(xi) |

Дn = max| yn+1(xi) - yn(xi) |)

n=1

д1 = 0,00932

n = 2



д2 = 2.615-10-4

Д1= 9.183*10-3

n=3

д3 = 2.688*10-5

Д2= 2.68*10-4

n=4



д4 = 5.405*10-7

Д3= 2.658*10-5

n = 5

д5 = 3.942*10-8

Д4= 5.411*10-7

Как показывают расчёты, с ростом n невязки стремятся к нулю.

Теперь обратимся к графикам:



Графики приближений слились, поэтому выполним их в искажённом масштабе (след. стр.)

*Графики выполнены в искажённом масштабе!

6.Результаты работы

В результате проделанной работы, мы получили следующие приближённые решения:

Первое:

yf1(x):=y1(x,a1)> 1,175201x - 0,267091x(1 - x)

Второе:

yf2(x):=y2(x,b1,b2)> 1,175201x -0.362747 x(1 - x)+0.191312x(1-x)2

Третье:

yf3(x):=y3(x,c1,c2,c3)> 1,175201x -0.367179x(1 - x)+ 0.213327x(1-x)2 -0.22015x(1-x)3

Четвёртое:

yf4(x):=y4(x,d1,d2,d3,d4)> 1,175201x -0.3678x(1 - x)+0.219442x(1-x)2 - 0.036272x(1-x)3+

+9.504139x10-3 x(1-x)4

Пятое:

yf5(x):=y5(x,e1,e2,e3,e4,e5)> 1,175201x +0.040634 x(1 - x)-0.019294x(1-x)2 - 1.012313x10-3

*x(1-x)3+2.408927*10-4 x(1-x)4+50356116*10-6x(1-x)5

В результате их сравнения были получены следующие невязки:

д1 = 9.32*10-3

д2 = 2.615-10-4

д3 =2.688*10-5

д4 = 5.405*10-7

д5 = 3.942*10-8

(отличия от точного решения)

Д1= 9.183*10-3

Д2= 2.68*10-4

Д3= 2.658*10-5

Д4= 5.411*10-7

(отличия приближений друг от друга)

Графики приближений и точного решения выглядят следующим образом:

*Графики выполнены в искажённом масштабе!

Выводы:

n	дi	Дi
1	0,00932
2	2.615-10-4	9.183*10-3
3	2.688*10-5	2.68*10-4
4	5.405*10-7	2.658*10-5
5	3.942*10-8	5.411*10-7

Таблица свидетельствует, что с ростом n невязки уменьшаются, т.е. приближения стремятся к точному решению задачи.

III. Приближённое решение краевой задачи для уравнения Эйлера методом конечных разностей

1. Описание метода конечных разностей

Для решения краевой задачи обыкновенного дифференциального уравнения (в нашем случае уравнения Эйлера) вида:

может быть использован метод конечных разностей. Суть его заключается в замене исходного дифференциального уравнения на соответствующее конечно - разностное. Сначала на отрезке [a,b] следует ввести сетку с шагом . Затем, на этой сетке определить сеточные функции pi = p(xi), qi = q(xi), fi = f(xi), отвечающие функциональным коэффициентам данного дифференциального уравнения. Считая y(x) точным решением данной краевой задачи, нужно обозначить через yi = y(xi) i - ую компоненту искомого каркаса приближённого решения. Фиксируя в исходном уравнении x = xi, с учётом наших обозначений следует перейти к равенствам:

уравнение эйлер функционал ритц

В каждом внутреннем узле сетки значения производных аппроксимируем конечноразностными отношениями:

Подставляем данные выражения в уравнение и получаем:

После приведения подобных, уравнение имеет вид:

В итоге, мы получили систему из (n-1) уравнений относительно yi. Используем краевые условия и решаем систему, решая, тем самым, исходную задачу в узлах yi.

2. Составление конечноразностного уравнения

Полученное нами уравнение Эйлера выглядит следующим образом:

Заданы, также, краевые условия y(0) = 0, y(1) = sh1.

Найдём конечноразностный аналог уравнения:

Данное равенство справедливо для k = 1..n-2 (где n = 10, 20, 100)

3. Описание метода прогонки

Пусть yk удовлетворяет уравнению akyk+1 - bkyk + ckyk-1 = fk (k = 1, …, n-1)

y0 = ч1 y1 + м1

yn = ч2 yn-1 + м2

Матрица этой системы будет иметь вид:

¦	1	- ч1	0	0	0	…	…	…	0	¦
¦	с1	-b1	a1	0	0	…	…	…	0	¦
¦	0	c2	-b2	a2	0	…	…	…	…	¦
¦	…	…	…	…	…	…	…	…	…	¦
¦	…	…	…	…	…	…	…	…	…	¦
¦	0	…	…	…	…	0	cn-1	-bn-1	an-1	¦
¦	0	…	…	…	…	0	0	-ч2	1	¦

a, b, c ? 0; ч или м могут равняться 0.

Система уравнений с такой матрицей (она называется тригональной) решается методом прогонки. Предположим, что решения для всех возможных значений целочисленного аргумента могут быть представлены в виде

yk = бk+1 y k+1 + в k+1 (k = 0…n-1) (1)

Отметим, что б и в не могут произвольными, даже если они обеспечивают это равенство. Дело в том, что элементы yk образуют столбец, являющийся решением краевой задачи для конечноразностного уравнения второго порядка. Поэтому для б и в существуют ограничения. Действительно, из (1) следует, что

yk-1 = бk y k + в k

akyk+1 - bkyk + ck(бk y k + в k) = fk

- bkyk + ckбk y k = fk - ck в k - akyk+1

Сравнивая, с предполагаемым решением:

(k = 0..n-1) (2)

(k = 0..n-1) (3)

В формуле (1) все предыдущие значения yk выражаются линейным образом через последующие, а в формуле (2) все последующие значения б выражаются через предыдущие. В формуле (3) все последующие значения в выражаются через предыдущие значения б и в.

Для первого краевого условия: б1 = ч1 (4)

Зная ч1 и используя формулу (2) можно найти все б. Аналогично в 1 = м1 и, с помощью формулы (3) можно найти все в.

Из второго краевого условия получаем:

¦ yn-1 = бn y n + в n

¦ yn = ч2 yn-1 + м2

yn = ч2 (бn y n + в n) + м2

(5)

В результате, алгоритм решения данной задачи методом прогонки сводится к следующему:

По (4) задаём б1 и в1, по (2) и (3) находим б и в (прогонка вперёд), по (5) находим yn, по формуле (1) находим все yk (прогонка назад)

4. Листинг программы на языке C++

#include <iostream.h>

#include <stdlib.h>

#include <iomanip.h>

#include <conio.h>

#include <math.h>

void main()

{ clrscr();

int n;

float x[110],toch[110],norm[110];

float nevyazka;

float b,c;

float h;

float y[110];

float alfa[110],beta[110];

cout<<"\ vvedite n->";

cin>>n;

cout<<"\n vvedite a->";

cin>>x[0];

cout<<"\n vvedite b->";

cin>>x[n];

h=(x[n]-x[0])/n;

cout<<h;

for(int i=0;i<n+2;i++) x[i]=(x[0]+i*h);

alfa[0]=0;

beta[0]=1;

b=2;

c=1-h*h;

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

{ alfa[i+1]=1/(b-c*alfa[i]);

beta[i+1]=c*beta[i]/(b-c*alfa[i]);

}

y[n]=sinh(1);

y[0]=1;

cout<<"\n yn";

for(int j=n-1;j>0;j--)

{y[j]=alfa[j+1]*y[j+1]+beta[j+1];

cout<<"\n"<<y[j];

}

cout<<"yt";

for(i=0;i<n+1;i++)

{toch[i]=sinh(x[i]);

cout<<"\n"<<toch[i]<<"\n" ;

norm[i]=fabs(y[i]-toch[i]);

}

if(n!=100)

{for(i=0;i<n+1;i++) cout<<"\n norm["<<i<<"] = "<<norm[i];

}

else for(int p=0;p<101;p+=10) cout<<"\n norm["<<p<<"] = "<<norm[p];

nevyazka=norm[0]; getch(); }

5. Расчёт приближений

n = 10

yn yт Дn

0.0000000000 0.0000000000 0.0000000000

0.1046919984 0.1001667500 0.0045252484

0.2093839967 0.2013360025 0.0080479942

0.3151229151 0.3045202934 0.0106026217

0.4229556734 0.4107523258 0.0122033476

0.5339396609 0.5210953055 0.0128443554

0.6491532051 0.6366535821 0.0124996230

0.7697061459 0.7585837018 0.0111224441

0.8967506188 0.8881059822 0.0086446366

1.0314921531 1.0265167257 0.0049754274

1.1752011936 1.1752011936 0.0000000000

n = 20

yn yт Дn

0.0000000000 0.0000000000 0.0000000000

0.1024224472 0.1001667500 0.0022556972

0.2053570067 0.2013360025 0.0040210042

0.3098298233 0.3045202934 0.0053095299

0.4168772906 0.4107523258 0.0061249648

0.5275563570 0.5210953055 0.0064610515

0.6429549861 0.6366535821 0.0063014040

0.7642028748 0.7585837018 0.0056191729

0.8924825343 0.8881059822 0.0043765521

1.0290408445 1.0265167257 0.0025241188

1.1752011936 1.1752011936 0.0000000000

n = 100

yn yт Дn

0.0000000000 0.0000000000 0.0000000000

0.1006174999 0.1001667500 0.0004507499

0.2021409800 0.2013360025 0.0008049775

0.3055851502 0.3045202934 0.0010648568

0.4119829204 0.4107523258 0.0012305946

0.5223957151 0.5210953055 0.0013004096

0.6379240630 0.6366535821 0.0012704809

0.7597185678 0.7585837018 0.0011348659

0.8889913678 0.8881059822 0.0008853856

1.0270282011 1.0265167257 0.0005114754

1.1752011936 1.1752011936 0.0000000000

yn - приближённое решение, yт - точное, Дn - невязки

6. Оценка невязок

Оценим невязки:

при n = 10 Дmax = 0.0128443554

при n = 20 Дmax = 0.0064610515

при n = 100 Дmax = 0.0013004096

Выводы

Так как невязки уменьшаются при увеличении числа точек на отрезке (измельчении шага), приближённое решение заданного дифференциального уравнения стремиться к точному.

Список литературы

1. Эльсгольц Л.Э. Дифференциальные уравнения и вариационное исчисление. - М., 1991.

2. Краснов М.Л., Макаренко Г.И., Киселёв А.И. Вариационное исчисление. - М., 1973.

3. Гельфанд И.М., Фомин С.В. Вариационное исчисление:М. 1961

4. Марчук Г.И. Методы вычислительной математики. - М.: Наука, 1977.

5. Березин И.С., Жидков Н.П. Методы вычислений. - М., 1987. т. II.

6. Самарский А.А., Гулин А.В. Численные методы: учебное пособие для вузов. - М.: Наука, 1989

7. Подбельский H.И. С++: Дрофа 2003

8. Дьяконов В.П. Энциклопедия MATHCAD 2001 и MATHCAD 11. С, 2004

Размещено на Allbest.ru