Основы информатики. Вычисление обыкновенных дифференциальных уравнений. Метод Рунге-Кутта (через for)

Для расчетов на профессиональном уровне наиболее часто используются средства из группы методов РУНГЕ-КУТТА. Для большинства задач наиболее оптимальным из этой группы является метод Рунге-Кутта четвертого порядка, который достаточно прост в реализации, имеет высокую точность и хорошую устойчивость. Мы по-прежнему рассматриваем дифференциальное уравнение первого порядка y'=f(x,y) с начальным условием y(x0)=y0. Для решения выбирается достаточно малый постоянный шаг изменения независимой переменной h так, что очередное значение xi есть xi-1+h=x0+ih, где i=1,2,3,.... Очередное значение искомой функции определяется из предыдущего по формуле:где коэффициенты k на каждом шаге определяются через значения функции f(x,y) при определенных значениях аргументов:Можно видеть, что на каждом шаге сначала вычисляются коэффициенты в той последовательности, в которой они указаны (поскольку они вычисляются один через другой), а затем определяется очередное значение функции.

Программный код:
procedure TForm1.Button3Click(Sender: TObject);
var
 x,y,h,k1,k2,k3,k4:real;
 i,n:integer;
begin
 n:=strtoint(edit1.Text);
 h:=1/n;
 x:=0;
 y:=1;
 paintbox1.Canvas.Pen.Color:=clgreen;
 paintbox1.Canvas.moveto(c.x+round(x*f.x),c.Y-10-round(y*f.y));
 for i:=1 to n do
  begin
  k1:=x*y*h;
  k2:=((x+h/2)*(y+k1/2))*h;
  k3:=((x+h/2)*(y+k2/2))*h;
  k4:=((x+h)*(y+k3))*h;
  x:=x+h;
  y:=y+(k1+k4+2*(k2+k3))/6;
  paintbox1.Canvas.lineto(c.x+round(x*f.x),c.Y-10-round(y*f.y));
  end;
end;

Если у кого то Delphi ниже 7-й версии, то удалите в исходнике слово XPMan в разделе Uses, и строчку XPManifest1: TXPManifest; в разделе type.

Скачать проект

Основы информатики. Вычисление обыкновенных дифференциальных уравнений.

Аппарат дифференциальных уравнений, несомненно, является наиболее мощным и наиболее распространенным средством описания процессов и явлений в самых разнообразных областях. Следовательно, именно дифференциальные уравнения являются наиболее часто используемыми инструментами математического моделирования. Именно поэтому дифференциальным уравнениям в данной главе уделено наибольшее внимание.
Предполагается, что читателю известны начальные сведения из теории обыкновенных дифференциальных уравнений (ОДУ). Освежив в памяти эти сведения, читатель вспомнит о том, что для ОДУ порядка выше первого (для определенности будем говорить об уравнениях второго порядка) возможны две принципиально различные постановки задачи решения уравнения. Если все начальные условия, определяющие, в конечном итоге, значения функции и ее производных, заданы в одной точке - на одном из концов интервала изменений независимой переменной, то говорят, что сформулирована задача Коши, или начальная задача. Если условия заданы на обоих концах интервала, на котором строится решение, то такая задача решения ОДУ называется краевой. Для уравнений первого порядка имеет смысл говорить только о начальной задаче, поскольку для них задается единственное условие. Основное наше внимание будет уделено решению задач Коши для ОДУ. Начнем с уравнения первого порядка, которое в общем виде можно представить в форме F(y',y,x)=0 c начальным условием y(x0)=y0. Далее будем считать, что уравнение может быть разрешено относительно производной так, что приводится к виду:

Все методы интегрирования ОДУ в задачах Коши сводятся к приближенному вычислению последующего значения yi в точке xi через предыдущее значение yi-1 в точке xi-1, при заданном из начального условия значении y0. В простейших случаях можно исходить из непосредственного определения понятия производной, переходя от бесконечно малых к конечным разностям:При этом нетрудно убедится, что из дифференциального уравнения следует:

Внимательный читатель заметит, что не мешало бы поставить индексы у x и y в аргументах функции f(y,x). Читатель, несомненно, прав. Вот вопрос о том, какие индексы поставить, весьма нетривиален. Вообще-то мы можем выбрать i или i-1 для независимой переменной и искомой функции произвольно. Все равно приближенное равенство должно быть справедливым. Обе возможности, действительно, можно реализовать. Если выбор за индексом i-1, то получаем формулу МЕТОДА ЭЙЛЕРА:

В противном случае имеем формулу НЕЯВНОЙ СХЕМЫ:

Задание:Решить дифференциальное уравнение y'=y на интервале от 0 до 2 с начальным условием y(0)=1 методом Эйлера и по неявной схеме. Сравнить оба численных решения с точным y=exp(x) путем построения графиков решений. Предусмотреть возможность изменения величины шага интегрирования, и убедиться в том, что с уменьшением шага точность улучшается.

Программный код для МЕТОДА ЭЙЛЕРА:
procedure TForm1.Button2Click(Sender: TObject);
var
x,y,k:real;
begin
paintbox1.Canvas.Pen.Color:=clred;
x:=0;
y:=1;
paintbox1.Canvas.MoveTo(c.X+round(x*f.x),c.y-round(y*f.y));
k:=strtofloat(edit1.Text);
while x<=2 do
begin
paintbox1.Canvas.LineTo(c.x+round(x*f.x),c.y-round(y*f.y));
x:=x+k;
y:=y+y*k;
end;
end;

Программный код для МЕТОДА НЕЯВНОЙ ФУНКЦИИ:
procedure TForm1.Button3Click(Sender: TObject);
var
x,y,k:real;
begin
paintbox1.Canvas.Pen.Color:=clgreen;
x:=0;
y:=1;
paintbox1.Canvas.MoveTo(c.X+round(x*f.x),c.y-round(y*f.y));
k:=strtofloat(edit1.Text);
while x<=2 do
begin
paintbox1.Canvas.LineTo(c.x+round(x*f.x),c.y-round(y*f.y));
x:=x+k;
y:=(y+y*k)*k+y;
end;
end;

Здесь c.x и c.y координаты начала трсчета(центр координатной системы), f.x и f.y фокусы масштабирования по соответствующим осям. Советую всем качать исходник. Если у кого то Delphi ниже 7-й версии, то удалите в исходнике слово XPMan в разделе Uses, и строчку XPManifest1: TXPManifest; в разделе type.

Скачать проект

Основы информатики. Вычисление определенных интегралов. Метод Симпсона.

Метод трапеций наиболее прост, но не является оптимальным по быстродействию. Ненамного более громоздким, но значительно более оперативным является МЕТОД СИМПСОНА. Не вдаваясь в детали, укажем, что в этом методе отдельные участки подынтегральной функции представляются не линейной, как в методе трапеций, а квадратичной интерполяцией. По этой причине метод Симпсона имеет второе название - метод парабол. В этой методике число разбиений области интегрирования n должно быть четным. Тогда формула Симпсона запишется в виде:Для вычисления интеграла с заданной точностью можно начать, например, с n=2 и далее действовать аналогично методу трапеций, удваивать число разбиений.
В данном разделе было бы естественным представить методы вычисления кратных (двойных, тройных и так далее) определенных интегралов. Однако, мы отложим этот вопрос на будущее. Дело в том, что применение к кратным интегралам концепций методов трапеций или парабол сталкивается с рядом серьезных проблем. Так, для двойных интегралов, область интегрирования следовало бы разбивать на прямоугольные или квадратные малые подобласти. При этом возникает достаточно нетривиальная задача правильного учета произвольной формы границы. Еще большие проблемы связаны с аппроксимацией "вырезаемого" подобластью участка телом с достаточно просто вычисляемым объемом. Хотя указанные проблемы в принципе решаемы, более эффективным оказывается применение здесь одного из приложений общей концепции методики статистических испытаний или метода Монте-Карло. Рассмотрению этой концепции посвящен последний раздел этой главы, где и описаны средства вычисления кратных интегралов.

Методом трапеций и методом Симпсона вычислить с заданной точностью определенный интеграл:

Программный код:

procedure TForm1.Button1Click(Sender: TObject);
var
p,x,h,a,b,xch,xnch:real;
n,i:integer;
begin
a:=0;
b:=pi;
n:=strtoint(edit1.Text);
h:=(b-a)/n;
x:=sin(a)+sin(b);
xch:=0;
xnch:=0;
for i:=1 to n do
if i mod 2=0 then
xch:=xch+sin(a+i*h)
else
xnch:=xnch+sin(a+i*h);
x:=(x+4*xnch+2*xch)*h/3;
edit2.Text:=floattostr(x);
p:=abs((2-x)/2);
edit3.Text:=floattostr(p);
end;

Скачать проект

Вся теория

Основы информатики. Вычисление определенных интегралов. Метод трапеций.

Необходимость в вычислении определенных интегралов с использованием численных методов возникает, в основном, в двух случаях. Во-первых, подынтегральная функция может быть задана таблично, например, как полученная в результате каких-либо измерений. Во-вторых, первообразная не может быть найдена в аналитическом виде, то есть, интеграл не является табличным. Как и все вычислительные методы, численное интегрирование является приближенным, и интеграл может быть получен с заданной точностью. Простейшим методом численного интегрирования является МЕТОД ТРАПЕЦИЙ. Речь идет о вычислении интегралаМетод применим, если пределы интегрирования конечны, а подынтегральная функция не имеет особенностей, хотя несложная предварительная процедура деления интервала интегрирования на несколько частей позволяет применить эту методику и к интегрированию функции, имеющей конечное число разрывов первого рода. В методе трапеций интервал [a,b] разбивается на n элементарных отрезков длиной

точками с координатой

На каждом из элементарных отрезков подынтегральная функция заменяется линейной функцией, а соответствующий участок площади под кривой y=f(x) заменяется трапецией так, как это показано на рис. 1.7.
Исходя из известной формулы для площади трапеции, нетрудно получить выражение для приближенного значения интеграла:
Эту формулу наиболее удобно использовать, если интегрируется функция, заданная таблично, причем, координаты xi могут быть расположены на интервале [a,b] произвольным образом (не обязательно равноудалены друг от друга).
Для вычисления интеграла от функции, заданной аналитически, целесообразно воспользоваться другим представлением формулы трапеций. В этом случае отрезок [a,b] может быть разбит на n равных отрезков длиной h. При этом предыдущая формула упрощается и принимает вид:

Методом трапеций и методом Симпсона вычислить с заданной точностью определенный интеграл:

Программный код:
procedure TForm1.Button1Click(Sender: TObject);
var p,x,h,a,b:real;
n,i:integer;
begin
a:=0;
b:=pi;
n:=strtoint(edit1.Text);
h:=(b-a)/n;
x:=(sin(a)+sin(b))/2;
for i:=1 to n do
x:=x+sin(a+i*h);
x:=x*h;
edit2.Text:=floattostr(x);
p:=abs((2-x)/2);
edit3.Text:=floattostr(p);
end;

Скчать проект

Вся теория

Основы информатики. Методы поиска корней уравнений. Метод хорд. Delphi.

Более быстродействующим, по сравнению с методом половинного деления, является способ поиска корня, названный МЕТОДОМ ХОРД. Графическая иллюстрация метода представлена на рисунке 1.6. Для построения алгоритма метода хорд необходима дополнительная информация о характере поведения функции на интервале поиска корня. Прежде всего, необходимо гарантировать, что на этом интервале функция монотонна (монотонно возрастает или монотонно убывает). Кроме того, на всем интервале не должен меняться характер выпуклости или вогнутости. Иными словами, на [a,b] не должны менять знак ни первая, ни вторая производные функции. Вообще говоря, даже и при нарушении этих условий метод хорд можно применять, но с использованием специальных приемов, на которых мы не будем останавливаться. Проще всего в сомнительных случаях просто сузить интервал до такого размера, на котором производные знак не меняют.

Рис. 1.6. Различные варианты в методе хорд

Из геометрии рисунка 1.6 (подобия треугольников) можно найти точку m пересечения хорды с осью абсцисс. Поскольку , то Дальнейшее построение алгоритма зависит от соотношения знаков первой и второй производных. Если знаки производных различны - левая часть рисунка 1.6, то новым правым концом интервала поиска корня становится точка m, то есть делается замена b на m. В противном случае, соответствующем двум вариантам правой части рисунка 1.6, делается замена на m. Итерационный процесс продолжается до достижения необходимой точности.
  Метод хорд действительно значительно более оперативен по отношению к методу деления отрезка пополам. В этом можно убедиться, реализовав оба метода для нахождения корня одного и того же уравнения на одном и том же отрезке и при одной и той же точности. Следует, однако, заметить, что выигрыш во времени будет иметь место только при оптимально написанной программе. В первую очередь, необходимо сократить до минимума количество вызовов функции f(x). Так в приведенном выше выражении для величины m по два раза фигурируют функции f(a) и f(b). Разумеется, функцию необходимо посчитать один раз, запомнить ее значение в какой-либо переменной и во второй раз использовать уже не вызов функции, а эту переменную.
  Рассмотренные методы не исчерпывают арсенала средств поиска корней уравнений. Существуют методы более оптимальные либо по времени, либо по расходованию памяти компьютера. Кроме того, имеются методики, оптимальные для конкретных видов функций, для которых решается уравнение. Так, например, существуют алгоритмы, позволяющие находить группы корней, если функция является алгебраическим полиномом (причем корни могут быть найдены в комплексной области). Алгоритмы эти выглядят, однако, довольно сложными и не могут быть рассмотрены здесь.

Решить уравнение sin(x)=0,5 на промежутке [0,1].

Программный код:
procedure TForm1.Button1Click(Sender: TObject);
var
 a,b,e,m,ms,s:real;
begin
  a:=0;
  b:=1;
  m:=1;
  e:=strtofloat(edit1.Text);
  s:=0;
  while abs((m-ms)/m)>e do
  begin
  ms:=m;
  m:=(a*(sin(b)-1/2)-b*(sin(a)-1/2))/((sin(b)-1/2)-(sin(a)-1/2));
  b:=m;
  s:=s+1;
  end;
  edit2.Text:=floattostr(m);
  edit3.Text:=floattostr(s*2);
end;

Скачать проект

Основы информатики. Методы поиска корней уравнений. Метод деления отрезка пополам. Delphi.

  В данном разделе речь идет о численных методах решения уравнения f(x)=0 в области действительных чисел. В общем, виде задача может быть сформулирована следующим образом. Необходимо найти такое значение x1, которое приближенно совпадает со значением x0, обращающим уравнение в тождество. Требуется, однако, существенно дополнить постановку задачи. Во-первых, поскольку уравнение может и не иметь решения, следует быть уверенным в том, что поставленная задача имеет смысл - нельзя поймать черную кошку в темной комнате, если ее там нет. Во-вторых, поскольку уравнение может иметь несколько корней, необходимо ограничить интервал поиска отрезком [a,b], для которого известно, что на нем имеется только один корень. Признаком этого является то, что f(a) и f(b) имеют разные знаки. Термин "приближенно совпадает" также требует пояснения. Все приближенные методы нахождения корня строятся на основе итераций, когда очередное значение так или иначе получается из предыдущего и все ближе подходит к точному значению корня. Теперь можно сформулировать условие прекращения этого итерационного процесса. Процесс поиска корня следует прекратить тогда, когда относительная разность двух последовательных приближений станет по абсолютной величине меньше заданной точности. Если максимальная допустимая ошибка задана величиной , а два последовательных приближения являются значениями xn-1 и xn, то условие прекращения итераций примет вид:

  Неприятность может произойти, если корнем уравнения является 0, однако это уже совсем особый случай, который рассматриваться не будет.
 
Рис. 1.5. Графическая иллюстрация метода половинного деления

  Наиболее наглядным, простым в реализации, однако, увы, не самым лучшим численным методом поиска корней, является МЕТОД ДЕЛЕНИЯ ОТРЕЗКА ПОПОЛАМ. Суть метода проще всего проиллюстрировать графически, с помощью рисунка 1.5.
  Прежде всего, находится точка , делящая отрезок [a,b], на котором ищется корень, пополам. Еще раз напомним, что корень на этом отрезке обязан быть. Далее, из этих двух половин выбирается та, на концах которой функция имеет разные знаки, то есть та, на которой имеется корень. Если эта половина является отрезком [c,b], как на нашем рисунке, то новым левым концом отрезка становится точка с. В противном случае точка с становится новым правым концом отрезка. Таким образом, после каждого такого выбора мы имеем новый отрезок [a,b], к которому снова применяем половинное деление. В этом методе наиболее наглядно выглядит условие выхода из цикла  А почему мы здесь обошлись без знака абсолютного значения? Автор надеется, что читатель сможет дать ответ на этот вопрос самостоятельно.

Решить уравнение sin(x)=0,5 на промежутке [0,1].

Программный код:

procedure TForm1.Button1Click(Sender: TObject);
var
 e,a,b,c,s:real;
begin
 a:=0;
 b:=1;
 s:=0;
 e:=strtofloat(edit1.Text);
 while abs((b-a)/b)>e do
  begin
  c:=(a+b)/2;
  if (sin(a)-1/2)*(sin(c)-1/2)>0 then
  a:=c
  else
  b:=c;
  s:=s+1;
  end;
 edit2.Text:=floattostr(c);
 edit3.Text:=floattostr(s*2);
end;

Скачать проект