Программная реализация решения обратной задачи методом наименьших квадратов

Министерство образования и науки Российской Федерации Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования ИРКУТСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ Институт недропользования Кафедра технологий геологической разведки ПОЯСНИТЕЛЬНАЯ ЗАПИСКА к курсовой работе по дисциплине Решение обратных задач в геофизике

Программная реализация решения обратной задачи методом наименьших квадратов Выполнил студент группы ГИС -10

Н.Б. Намдаков Нормоконтроль Е. В. Агеенков Иркутск 2014 г.

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

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

Для того, чтобы получить представления о свойствах объекта, необходимо создать его модель.

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

Процесс перехода от модели к полю называют решением прямой задачи. Переход от значения поля к параметрам модели среды — решением обратной задачи. Одним из простейших вариантов решения обратной задачи является подбор такой модели, которая дала бы теоретическое поле, совпадающее или близкое к наблюдаемому.

Прямая задача, как правило, имеет единственное решение. Заданной модели при заданных условиях наблюдения соответствует единственное поле. В обратной задаче — одному и тому же полю может соответствовать множество моделей. Поэтому при решении задач возникает вопрос: какой ответ мы получили единственный или один их множества и какой из множества ответов наиболее близок к реальному. Прямые задачи являются, как правило, устойчивыми. Обратные задачи очень часто оказываются неустойчивыми, т. е. небольшие искажения в данных наблюдений могут приводить к значительным погрешностям в параметрах модели.

1. Теоретическая часть

1.1 Программные средства

Наиболее удобной, в смысле программной реализации, средой программирования является Delphi. Среда визуального программирования Delphi, разработанный компанией Borland. Среда Delphi включает в себя полный набор визуальных инструментов для скоростной разработки приложений, поддерживающий разработку пользовательского интерфейса и подключение к корпоративным базам данных. VCL — библиотека визуальных компонент, включает в себя стандартные объекты построения пользовательского интерфейса, объекты управления данными, графические объекты, объекты мультимедиа, диалоги и объекты управления файлами, управление DDE и OLE.

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

1.2 Физическая модель

Апроксимация функций методом наименьших квадратов

При решении многих практических задач часто приходится вычислять значения каких-то функциональных зависимостей y = f (x).

При этом, как правило, имеют преобладающее место две ситуации.

1. Явная зависимость между х и y на [a, b] отсутствует, а имеется только таблица экспериментальных данных {x_i, y_i}, и возникает необходимость определения y = f (x) на интервале [x_i, x_i/2] [a, b]. К этой задаче относится также уточнение таблиц экспериментальных данных.

2. Зависимость y = f (x) известна и непрерывна, но настолько сложна, что не пригодна для практических расчетов. Стоит задача упрощения вычисления значений y = f (x) и ее характеристик (и т.д.). Поэтому, с точки зрения экономии времени и материальных ресурсов, приходят к необходимости построения какой-то другой функциональной зависимости y = F (x), которая была бы близка к f (x) по основным ее параметрам, но более проста и удобна в реализации при последующих расчетах, т. е. ставится задача о приближении (аппроксимации) в области определения y = f (x). Функцию y = F (x) называют аппроксимирующей.

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

Метод наименьших квадратов

Пусть таблично задана функция. Определить аппроксимационный многочлен вида:

(1)

Для определения коэффициентов (где i=1,2…m) используют аппроксимацию методом наименьших квадратов. Функционал находится по формуле

где n — количество пар значений аргумента x_i и функции y_i. Т.к. S должен быть минимален, то первые частные производные от S по всем коэффициентам полинома должны равняться нулю.

Вычислим их и приравняем к нулю:

(1)

После преобразования система слегка упростится:

(2)

Если полином первой степени, то 2 уравнения, если шестой степени, то 7 уравнений. Введём следующие обозначения:

С учётом этих обозначений система (2) перепишется следующим образом:

(3)

Решив эту систему одним из известных методов определится ряд значений, с помощью которого и строится многочлен (1).

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

1.3 Алгоритм решения

аппроксимация квадрат программирование полином

1. Вводим исходные данные Xi и Yi из файла.

2.Вычисляем коэффициенты Т [ссылка на теорию по методу], где

3. Вычисляем коэффициенты С, где

4. Формируем расширенную матрицу системы уравнений, т. е. из системы:

с учетом введённых обозначений получаем систему:

5.Приводим расширенную матрицу системы уравнений к треугольному виду (прямой ход метода Гаусса).

6.Вычисляем коэффициенты полинома (обратный ход метода Гаусса).

7.Выводим коэффициенты полинома (в зависимости от степени полинома используем массивы соответствующей размерности и присваиваем элементам значения, полученные после решения системы методом Гаусса).

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

1.4 Блок-схема

2. Практическая часть

2.1 Решение задачи в Delphi

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

На форме у нас 3 компоненты: TMemo, и 4 компоненты TButton (рис. 1).

Рис. 1. Главная форма

Заходим в инспектор объектов компонента TButton1, на свойстве Caption пишем: «Загрузить данные». Таким же образом меняем свойство Caption на «Выполнить», «График» и «Выход» соответственно.

Создаем вторую форму. Кидаем на форму компоненты: TPanel, TChart и TButton. Главная рабочая компонента это TChart, в основном работаем с ней. Щелкните два раза, откроется следующее диалоговое окно (рис. 2).

Рис. 2. Компонента TChart

Далее нажимаем на кнопку «Add…» создаем новые линии. В данной задаче у нас 4 линии: полиномы 3-ей, 5-ой и 7-ой степени, а также исходные данные. Исходные данные необходимы для сравнения, сильно ли исказился наш график.

Дальше пишем процедуру программы [Приложения].

2.2 Графическое представление результатов

Нахождение коэффициентов аппроксимирующих полиномов 3-ей, 5-ой и 7-ой степени (рис. 3).

Рис. 3. Выполнение программы

Графическое представление нахождения коэффициентов аппроксимирующих полиномов (рис. 4).

Рис. 4. График функции полиномов и исходных данных

Заключение

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

С помощью данного алгоритма можно вычислить коэффициенты аппроксимирующего полинома и другой степени (не меньшей 1 и не большей 8) при другом количестве экспериментальных точек.

Использованная литература

1. Delphi для «чайников». Нейл Дж. Рубенкинг. Киев — Москва: Диалектика, 1997.

2. http://www.delphi-manual.ru/

3. Метод наименьших квадратов (МНК) 2010 г. — (http://www.cleverstudents.ru/articles/mnk.html)

4. Аппроксимация функции полиномом методом наименьших квадратов — (http://alexeypetrov.narod.ru/C/sqr_less_about.html)

Приложение

unit Unit1;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, StdCtrls, Series, TeEngine, TeeProcs, Chart, ExtCtrls, buttons, Math, Unit2;

type

TForm1 = class (TForm)

Memo1: TMemo;

Memo2: TMemo;

Button1: TButton;

Memo3: TMemo;

Button2: TButton;

Button3: TButton;

Button4: TButton;

Label1: TLabel;

Label2: TLabel;

procedure Button1Click (Sender: TObject);

procedure Button2Click (Sender: TObject);

procedure Button3Click (Sender: TObject);

procedure FormCreate (Sender: TObject);

procedure Button4Click (Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;

var

Form1: TForm1;

p3: array [0. 8] of extended;

p5: array [0. 8] of extended;

p7: array [0. 8] of extended;

x: array [0. 39] of extended;

f: array [0. 39] of extended;

implementation

{$R *.dfm}

procedure Polynom (n, m: integer);

var

T: array [0. 16] of extended;

C, A: array [0. 8] of extended;

B: array [0. 8, 0. 9] of extended;

DataX, DataY: text;

StrXi, StrYi: string;

i, ii, j, jj, k, s, Code: integer;

Xi, Yi, Y, Bik, Delta: extended;

begin

ii := 0;

for i := 0 to 16 do

T[i] := 0;

for i := 0 to 8 do

begin

C[i] := 0;

A[i] := 0

end;

for i := 0 to 8 do

for j := 0 to 9 do

B[i, j] := 0;

assign (DataX, 'Xi.txt');

assign (DataY, 'Yi.txt');

reset (DataX);

reset (DataY);

for i := 1 to m do

begin

readln (DataX, StrXi);

readln (DataY, StrYi);

val (StrXi, Xi, Code);

val (StrYi, Yi, Code);

x[ii] := Xi;

f[ii]: = Yi;

ii := ii + 1;

for j := 1 to 2 * n do

T[j] := T[j] + exp (j * ln (Xi));

for j := 0 to n do

C[j] := C[j] + Yi * exp (j * ln (Xi));

end;

T[0] := m;

close (DataX);

close (DataY);

for i := 0 to n do

for j := 0 to n do

B[i, j] := T[j + i];

for i := 0 to n do

B[i, n + 1] := C[i];

for k := 0 to n — 1 do

for i := k to n do

begin

Bik := B[i, k];

for j := k to n + 1 do

if i = k then

B[i, j] := B[i, j] / Bik

else

B[i, j] := B[i, j] / Bik — B[k, j];

end;

Form1.Memo3.Lines.Add ('Коэффициенты полинома степени: ' + IntToStr (n));

for i := n downto 0 do

A[i] := (B[i, n + 1] - B[i, 1] * A[1] - B[i, 2] * A[2] - B[i, 3] * A[3] - B

[i, 4] * A[4] - B[i, 5] * A[5] - B[i, 6] * A[6] - B[i, 7] * A[7] - B

[i, 8] * A[8]) / B[i, i];

for i := 0 to n do

begin

if (n = 3) then

p3[i] := A[i]

else if (n = 5) then

p5[i] := A[i]

else

p7[i] := A[i];

Form1.Memo3.Lines.Add ('A[' + IntToStr (i) + ']=' + FloatToStr (A[i]));

end;

Form1.Memo3.Lines.Add ('Среднеквадратичное отклонение:');

reset (DataX);

reset (DataY);

Delta := 0;

for i := 1 to m do

begin

readln (DataX, StrXi);

readln (DataY, StrYi);

val (StrXi, Xi, Code);

val (StrYi, Yi, Code);

Y := 0;

for j := 0 to n do

Y := Y + A[j] * exp (j * ln (Xi));

Delta := Delta + sqr (Y — Yi);

end;

Delta := sqrt (Delta / m);

Form1.Memo3.Lines.Add ('Дельта = ' + FloatToStr (Delta));

close (DataX);

close (DataY);

end;

procedure TForm1. Button1Click (Sender: TObject);

begin

Button2.Enabled := true;

Memo1.Lines.LoadFromFile ('Xi.txt');

Memo2.Lines.LoadFromFile ('Yi.txt');

end;

procedure TForm1. Button2Click (Sender: TObject);

begin

Polynom (3, 40);

Polynom (5, 40);

Polynom (7, 40);

Button3.Enabled := true;

end;

procedure TForm1. Button3Click (Sender: TObject);

var i: integer;

begin

Form2.Show;

for i := 0 to 39 do

begin

Form2.Series1.AddXY (x[i], (p3[0] + p3[1] * x[i] + p3[2] * sqr (x[i])

+ p3[3] * power (x[i], 3)));

Form2.Series2.AddXY (x[i], (p5[0] + p5[1] * x[i] + p5[2] * sqr (x[i])

+ p5[3] * power (x[i], 3) + p5[4] * power (x[i], 4)

+ p5[5] * power (x[i], 5)));

Form2.Series3.AddXY (x[i], (p7[0] + p7[1] * x[i] + p7[2] * sqr (x[i])

+ p7[3] * power (x[i], 3) + p7[4] * power (x[i], 4)

+ p7[5] * power (x[i], 5) + p7[6] * power (x[i], 6)

+ p7[7] * power (x[i], 7)));

Form2.Series4.AddXY (x[i], f[i]);

end;

end;

procedure TForm1. FormCreate (Sender: TObject);

begin

Memo1.Clear;

Memo2.Clear;

Memo3.Clear;

Button2.Enabled := false;

Button3.Enabled := false;

end;

procedure TForm1. Button4Click (Sender: TObject);

begin

close;

end;

end.

unit Unit2;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, StdCtrls, Series, TeEngine, TeeProcs, Chart, ExtCtrls, buttons, Math;

type

TForm2 = class (TForm)

Panel1: TPanel;

Chart1: TChart;

Series2: TLineSeries;

Series3: TLineSeries;

Series4: TPointSeries;

Series1: TLineSeries;

Button1: TButton;

procedure Button1Click (Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;

var

Form2: TForm2;

implementation

{$R *.dfm}

procedure TForm2. Button1Click (Sender: TObject);

begin

Form2.Close;

end;

end.