Наши проекты:
Журнал · Discuz!ML · Wiki · DRKB · Помощь проекту |
||
ПРАВИЛА | FAQ | Помощь | Поиск | Участники | Календарь | Избранное | RSS |
[3.137.218.230] |
|
Сообщ.
#1
,
|
|
|
Решение уравнений
{Автор: Ozzя} var a,b,x,e:real; { Нелинейное уравнение (x-2)*lg(x+11)-1 } function f(x:real):real; begin f:=(x-2)*lg(x+11)-1; end; { Функция, определяющая знак числа } function sgn(x:real):integer; begin sgn:=0; if x<0.0 then sgn:=-1; if x>0.0 then sgn:=1; end; { Процедура, реализующая метод половинного деления } procedure dich (var a,b,x,e:real); var i:integer; r:real; begin i:=sgn(f(a)); { Повторяем до тех пор, пока интервал [a,b] не станет меньше } { заданной погрешности е } while b-a>e do begin { Определяем середину отрезка } x:=(a+b)/2; { Далее делаем выбор, какую из частей отрезка взять } { для дальнейшего уточнения корня } { Если левая часть уравнения f(x) есть непрерывная функция } { аргумента х, то корень будет находиться в той половине } { отрезка, на концах которой f(x) имеет разные знаки. } r:=f(x); if sgn(r)=i then a:=x else b:=x end end; begin write('Введите значения интервала A,B? '); readln(a,b); write('Введите требуемое значение точности E? '); readln(e); dich(a,b,x,e); writeln('X= ',x); end. Это сообщение было перенесено сюда или объединено из темы "Решение уравнений" Это сообщение было перенесено сюда или объединено из темы "Заготовка для "Численных методов"" |
Сообщ.
#2
,
|
|
|
Численное решение уравнений итерационными методами.
Автор: DonNTU. e-moe aka Labinskiy Nikolay © 2005 Данная статья посвящена нахождению корней уравнения F(x)=0, (1) где функция F(x) может быть алгебраической либо трансцендентной и должна удовлетворять условию дифференцируемости. Как правило, численное решение уравнений состоит из двух этапов: нахождение приближенного значения корня и его уточнение до заданной точности. Начальное приближение часто известно из физических соображений либо находится специальными методами, например, графически. Подробно будет рассмотрен второй этап решения уравнений: нахождение корня с заданной точностью различными итерационными методами. при этом условия сходимости принимают вид: Применим метод Ньютона-Рафсона согласно формуле (10), при этом вычисление F(x) будем осуществлять по правилу Горнера с использованием рекуррентных формул: b[m] = a[m] b[j] = a[j] + x[n]*b[j+1], j=m-1, ..., 0 (11) Таким образом находим F(x) = b[0]. F'(x) представляет собой многочлен степени m-1. Воспользовавшись для его вычисления теми же рекуррентными формулами, имеем: c[m] = b[m] c[j] = b[j] + x[n]*c[j+1], j=m-1, ..., 1 (12) и соответственно F'(x) = c[1]. Подставляя найденные значения F(x) и F'(x) в формулу (10) для метода Ньютона-Рафсона, получаем x[n+1] = x[n] - b[0]/c[1], (13) где b[o] и c[1] вычислены по формулам (11) и (12). Тихо и незаметно мы пришли к методу Бриге-Виетта { DonNTU. e-moe aka Labinskiy Nikolay (c) 2005 visit www.sources.ru } {$N+,G+,X+} const Eps = 0.00001; { Точность вычислений } anMax = 3; { Степень полинома } { Коэффициенты полинома F(x) = a[0]+a[1]*x+a[2]*x^2+ ... +a[n]*x^n } An: array[0..anMax] of double = (-5.372,1.2493,0.559,-0.13); { a[0] a[1] a[2] a[3] } type { Описание функции для передачи ее в кач-ве параметра функциям } TFunc = function(x: double; px: byte; d: shortint): double; var x0,x1,x2: double; { Переменные для вычисления корней } function Func(x: double; px: byte; d: shortint): double; far; (* Функция для вычисления полинома в точке 'x' x - в этой точке будет происходить расчет px - если 1, то результат будет = F(x)+x d - если 0, то F(x), 1 - F'(x), -1 - F(x)/F'(x) *) var bn, { Это будет F(x) } cn: double; { Это будет F'(x) } j: byte; { Счетчик цикла ;) } begin if px=1 then { Если px=1, то добавляем к рез-ту x } An[1]:=An[1]+1; { Далее идет расчет по правилу Горнера } bn:=An[anMax]; cn:=bn; for j:= anMax-1 downto 1 do begin bn:=An[j]+x*bn; cn:=bn+x*cn; end; bn:=An[0]+x*bn; { В зависимости от d возвращаем результат } case d of -1: Func:=bn/cn; 0: Func:=bn; 1: Func:=cn; else Func:=bn; end; { Забираем обратно свой x } if px=1 then An[1]:=An[1]-1; end; function Sign(x: double): shortint; { Функция определения знака } begin if x < - Eps then Sign:=-1 else if x > Eps then Sign:=1 else Sign:=0; end; function Get_X0(f: TFunc): double; (* Поиск начального приближения: Идет начальный просмотр функции и определение двух точек, между которыми функция меняят знак (проходит через 0). Где-то между этими точками лежит корень. *) const LookFrom = -10.0; { Поиск начинается с этой точки } step = 1; { Шаг, с которым производится поиск } begin x0:=LookFrom; x1:=x0+eps; while Sign(f(x1,0,0)) = Sign(f(x0,0,0)) do begin x0:=x1; x1:=x1+step; end; Get_X0:= (x1+x0)/2; end; function Metod1(f: TFunc; x0: double): double; var cnt: word; function TestFunc(x: double): boolean; begin if abs(f(x,1,1)) >= 1 then begin writeln(' Ошибка: Невозможно расчитать результат!'); writeln(' Модуль производной в точке ',x:0:6,' больше 1',#10); TestFunc:= false; end else TestFunc:= true; end; begin Writeln('Метод последовательных приближений:'); if TestFunc(x0) then x1:=f(x0,1,0) else exit; { x2:=f(x1,0)+x1;} cnt:=1; while abs(x1-x0) > eps do begin x0:=x1; if TestFunc(x0) then x1:=f(x0,1,0) else exit; inc(cnt) end; Metod1:=x1; Writeln('x=', x1:0:6); Writeln('Метод сошелся на ',cnt,' шаге.',#10); end; function Metod2(f: TFunc; x0: double): double; var alph: double; cnt: word; begin Writeln('Усовершенствованный метод последовательных приближений:'); x1:=f(x0,1,0); alph:=1/(1-f(x1,0,0)/(x1-x0)); x2:=x1+alph*f(x1,0,0); cnt:=2; while abs(x2-x1) > eps do begin x0:=x1; x1:=x2; alph:=1/(1-f(x1,0,0)/(x1-x0)); x2:=x1+alph*f(x1,0,0); inc(cnt) end; Metod2:=x2; Writeln('x=', x2:0:6); Writeln('Метод сошелся на ',cnt,' шаге.',#10); end; function Metod3(f: TFunc; x0: double): double; var cnt: word; begin Writeln('Метод Бирге-Виетта:'); x1:=x0-f(x0,0,-1); cnt:=1; while abs(x1-x0) > Eps do begin x0:=x1; x1:=x0-f(x0,0,-1); inc(cnt); end; Metod3:=x1; Writeln('x=', x1:0:6); Writeln('Метод сошелся на ',cnt,' шаге.',#10); end; begin x0:=Get_X0(Func); Writeln; Writeln('F(x)= -5.372 + 1.2493*x + 0.559*x^2 - 0.13*x^3'); Writeln('Eps = ',eps:0:5,' x0= ',x0:0:6,#10); Metod1(Func,x0); Metod2(Func,x0); Metod3(Func,x0); Write('Press Enter to exit'); Readln; end. Это сообщение было перенесено сюда или объединено из темы "Решение уравнений" |
Сообщ.
#3
,
|
|
|
Численное решение систем линейных алгебраических уравнений Данная статья посвящена нахождению корней систем линейных алгебраических уравнений. Методы численного решения таких систем подразделяются на два типа: прямые (конечные) и итерационные (бесконечные). Оба типа полезны и удобны для практических вычислений и каждый из них имеет свои преимущества и недостатки. 1. Метод исключений (метод Гаусса). Данный метод является наиболее известным и широко применяемым. Рассмотрим систему из n уравнений с n неизвестными. Обозначим неизвестные через x[1], x[2], ... , x[n] и запишем систему в следующем виде: a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1] a[2,1]*x[1] + a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2] .................................................... a[n,1]*x[1] + a[n,2]*x[2] + ... + a[n,n]*x[n] = b[n] Предполагается, что в силу расположения уравнений, a[i,i] <> 0. Если это не так, то, меняя уравнения местами, добиваемся выполнения этого условия. Введем n-1 множителей: m[i] = a[i,1]/a[1,1] , где i = 2, 3, ... , n и вычтем из каждого i-го уравнения первое, умноженное на m, обозначая a'[i,j] = a[i,j] - m[i]*a[1,j] b'[i] = - m[i]*b[i] i = 2, 3, ... , n; j = 1, 2, ... , n. Для всех уравнений, начиная со второго, получим: a[i,1] = 0 , где i = 2, 3, ... , n. Преобразованная система уравнений запишется в следующем виде: a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1] 0 + a'[2,2]*x[2] + ... + a'[2,n]*x[n] = b'[2] .................................................... 0 + a'[n,2]*x[2] + ... + a'[n,n]*x[n] = b'[n] Продолжая таким же образом, исключаем x[2] из последних n-2 уравнений, затем x[3] из последних n-3 и т.д. На некотором k-м этапе мы исключаем x[k] с помощью множителей: m^(k-1)_[i] = a^(k-1)_[i,k]/a^(k-1)_[k,k] , i = k+1, ... , n, не забывая, что a(k-1)[k,k] <> 0. Примечание, ^(k-1) обозначает на k-1й итерации _[i,k] - индекс эл-та. Тогда a^(k)_[i,j] = a^(k-1)_[i,j] - m^(k-1)_[i]*a^(k-1)_[k,j] b^(k)_[i] = - m^(k-1)_[i]*b^(k-1)_[k] i = k+1, k+2, ... , n; j = k, k+1, ... , n. Окончательно, треугольная система уравнений записывается: a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1] a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2] .................................................... a[n,n]*x[n] = b[n] Обратная подстановка для нахождения значений неизвестных задается формулами: x[n] = b^(n-1)_[n] / a^(n-1)_[n,n] x[n-1] = (b^(n-2)_[n-1] - a^(n-2)_[n-1,n]*x[n]) / a^(n-2)_[n-1,n-1] ................................................................... x[j] = ( b^(j-1)_[j] - a^(j-1)_[j,n]*x[n] - ... - - a^(j-1)_[j,j+1]*x[j+1] ) / a^(j-1)_[j,j] j = n-2, n-3, ... , 1. Перед началом процесса исключения очередного неизвестного может понадобиться переставить уравнения в системе, чтобы избежать деления на нуль. Кроме этого, при перестановке уравнений необходимо добиваться того, чтобы | a^(k-1)_[k,k] | >= | a^(k-1)_[i,k] | Поэтому перед началом процесса исключения очередного неизвестного необходимо переставить уравнения так, чтобы максимальный по модулю коэффициент при x[k] попал на главную диагональ (данный способ решения часто называется методом главного элемента). 2. Итерационный метод Гаусса-Зейделя. Данный метод отличается простотой и легкостью программирования, малой ошибкой округления, но сходится при выполнении некоторых условий. Рассмотрим все ту же систему из n уравнений с n неизвестными. По прежнему полагаем, что a[i,i] <>0 для всех i. Тогда k-е приближение к решению будет задаваться формулой: x^(k)_[i] = ( b[i] - a[i,1]*x^(k)_[1] - ... - a[i,j-1]*x^(k)_[i-1] - - a[i,j+1]*x^(k-1)_[i+1] - a[i,n]*x^(k-1)_[n] ) / a[i,i] i = 1, 2, ... , n. Как известно, любой итерационный процесс, как и настоящий мужик, должен имееть свой конец Поэтому вычисления нужно продолжать до тех пор, все x^(k)_[i] не станут достаточно близки к x(k-1)[i], т.е. пока для всех i не будет выполнятся неравенство max| x^(k)_[i] - x^(k-1)_[i] | < Eps, где Eps некоторое положительное число, задающее точность вычислений. Итерационный метод Гаусса-Зейделя сходится, если удовлетворяется достаточное условие сходимости: для всех i модуль диагонального коэффициента должен быть не меньше суммы модулей остальных коэффициентов данной строки | a[i,i] | >= | a[i,1] | + | a[i,2] | + ... + | a[i,j-1] | + + | a[i,j+1] | + ... + | a[i,n] | , а хотя бы для одной строки это неравенство должно выполнятся строго | a[i,i] | > | a[i,1] | + | a[i,2] | + ... + | a[i,j-1] | + + | a[i,j+1] | + ... + | a[i,n] |. 3. Сравнение методов. Метод исключений имеет то преимущество, что он конечен и теоретически, с его помощью можно решить любую невырожденную систему. Итерационный метод сходится только для специальных систем уравнений, однако когда итерационные методы сходятся, они, обычно, предпочтительнее: А вот и исходник по теме: { Writeln('Labinskiy Nikolay aka e-moe (c) 2005'); Writeln('See you later on forum.sources.ru ;)'); } {$N+,G+,X+} const Comments: boolean = false; { нужно ли печатать пошаговые рез-ты расчета } Eps = 0.00001; { Точность вычислений } n_max = 4; { Кол-во ур-й и неизвестных } const_AnB: array[1..n_max,1..n_max+1] of double = (( -8, 2, 17, -5, -119.97), { -8*X[1] + 2*X[2] + 17*X[3] - 5*X[4] = -119.97 и т.д. ... } ( 4, -22, 6, 5, 55.73), ( 15, 3, -5, -5, 18.77), ( -4, -4, 5, 14, -79.42)); { Описание типов матрицы уравнеиний и массива найденных Х } Type TMatr = array[1..n_max,1..n_max+1] of double; Type TXMarr = array[1..n_max] of double; procedure PrintX(const X: TXMarr; n: byte; cnt: word); { Печатает найденные Х } var k, { Позиция заданного корня в массиве } i: byte; { Номер текущего корня для вывода } begin if n = 0 then exit; k:=low(X); if cnt = 0 then begin write(' N x1'); for i:=2 to n do write(' x',i); writeln; end; write(cnt:2); for i:=1 to n do begin Write(' ',X[k]:3:5); inc(k); end; writeln; { При решении систем, с большим кол-вом неизвестных, возможно, придется переделать эту процедуру для корректного вывода всех решений. } end; { PrintX } procedure PrintMatr(var AnB: TMatr; n: byte); { Печатает заданную матрицу } var i,j: byte; begin for i:=1 to n do begin for j:=1 to n+1 do write(AnB[i,j]:3:5,' '); writeln; end; writeln; { При решении систем, с большим кол-вом неизвестных, возможно, придется переделать эту процедуру для корректного вывода всех решений. } end; { PrintMatr } function Normalize(var AnB: TMatr; n: byte): byte; { Располагает на главной диагонали наибольшие элементы в столбцах Если один из них равен нулю,то систему решить не получится ... Возвращает 0 если систему заданными методами решить не получится, 1 - с-му можно решить только методом Гаусса 2 - с-му можно решать любым методом } var max: double; { Переменная для поиска макс. эл-та } imax: byte; { Переменная для поиска макс. эл-та } tmp: double; { Переменная для обмена элементов } i,j: byte; { Счетчики циклов } begin Normalize:=2; for j:=1 to n do begin max:=Abs(AnB[1,j]); imax:=1; for i:=2 to n do if Abs(AnB[i,j]) > max then begin max:=Abs(AnB[i,j]); imax:=i; end; if max < Eps then begin Normalize:=0; Writeln('Error: a[i,i]=0!'); exit; end else if j <> imax then for i:=1 to n+1 do begin tmp:=AnB[j,i]; AnB[j,i]:=AnB[imax,i]; AnB[imax,i]:=tmp; end; end; for i:=1 to n do begin tmp:=0; for j:=1 to n do if i <> j then tmp:=tmp+Abs(AnB[i,j]); if Abs(AnB[i,i]) < tmp then Normalize:=1; end; end; { Normalize } procedure Metod1(const Matr; n: byte); var k : byte; M_ : TMatr absolute Matr; AnB : TMatr; X,M : TXMarr; function Calc(k: byte): boolean; var i,j : byte; begin Calc:= true; for i:=k+1 to n do begin M[i]:= AnB[i,k]/AnB[k,k]; if M[i] > 1 then begin Calc:= false; exit; end; end; for i:=k+1 to n do for j:=k to n+1 do AnB[i,j]:=AnB[i,j]-M[i]*AnB[k,j]; if Comments then PrintMatr(AnB,n); end; { Calc } function Summ(j: byte): double; var i : byte; res : double; begin res:=AnB[j,n+1]; for i:=n downto j+1 do res:=res - AnB[j,i]*X[i]; Summ:=res; end; { Summ } begin Writeln('Метод Гаусса:'); AnB:=M_; if Normalize(AnB,n) = 0 then begin writeln('Error: Невозможно решить систему уравнений!'); exit; end; for k:=1 to n-1 do if not Calc(k) then begin Writeln('Error: Невыполняется условие | a^(k-1)_[k,k] | >= | a^(k-1)_[i,k] | !'); exit; end; X[n]:=AnB[n,n+1]/AnB[n,n]; for k:=n-1 downto 1 do X[k]:=Summ(k)/AnB[k,k]; PrintX(X,n,0); Writeln; end; { Metod1 } procedure Metod2(const Matr; n: byte); var M_ : TMatr absolute Matr; AnB : TMatr; X : TXMarr; i,j: byte; tmp: double; delta: double; cnt: word; function Summ(i: byte): double; var j: byte; res: double; begin res:=AnB[i,n+1]; for j:=1 to n do if i <> j then res:=res-AnB[i,j]*X[j]; Summ:=res; end; { Summ } begin writeln('Метод Гаусса-Зейделя:'); AnB:=M_; For i:=1 to n do X[i]:=0; if Normalize(AnB,n) <> 2 then begin writeln('Error: !'); exit; end; delta:=1; cnt:=0; while delta > Eps do begin if Comments then PrintX(X,n,cnt); delta:=0; for i:=1 to n do begin tmp:=Summ(i)/AnB[i,i]; if delta < Abs(tmp-X[i]) then delta:=Abs(tmp-X[i]); X[i]:=tmp; end; inc(cnt); end; PrintX(X,n,cnt); writeln('Заданная точность вычислений достигнута на ',cnt,' шаге.'); Writeln; end; { Metod2 } begin Writeln(#10,'Численное решение систем линейных алгебраических уравнений:',#10); Metod1(const_AnB,n_max); Writeln('Press Enter to continue'); readln; Metod2(const_AnB,n_max); Writeln('Press Enter if you wanna DOS ;)'); readln; end. Это сообщение было перенесено сюда или объединено из темы "Решение уравнений" |
Сообщ.
#4
,
|
|
|
Методы решения дифференциальных уравнений.
Данная статья посвящена численным методам решения дифференциальных уравнений. Мы будем рассматривать методы решения одного обыкновенного дифференциального уравнения первого порядка с одним начальным условием: y' = f(x,y) (1) y(x[0]) = y[0] (2) Методы, которые здесь рассматриваются, легко обобщаются для системы уравнений первого порядка. А уравнения высших порядков можно свести к системе уравнений первого порядка. В основном существуют два широких класса методов: Как и во многих других случаях, эти два класса методов придется сочетать разумным образом, учитывая их достоинства и недостатки. 1. Методы Рунге-Кутта. Методы Рунге-Кутта обладают следующими отличительными свойствами: Именно благодаря 3) эти методы удобны для практических вычислений, однако для вычисления одной последующей точки решения нам придется вычислять f(x,y) несколько раз при различных x и y. 1.1. Метод Эйлера. Этот метод, один из самых старых и широко известных, описывается формулой: y[m+1] = y[m] + h*f(x[m],y[m]). (3) Найденное по формуле (3) решение согласуется с разложением в ряд Тейлора вплоть до членов порядка h, т.е. данный метод является методом Рунге-Кутта первого порядка. Этот метод имеет довольно большую ошибку приближения, кроме того, он очень часто оказывается неустойчивым - изначально малая ошибка (происходящая от приближения, округления или исходных данных) увеличивается с ростом x. Для вычисления y[m+1] метод Эйлера использует наклон касательной только в точке x[m],y[m]. Этот метод можно усовершенствовать множеством различных способов. Рассмотрим два из них. 1.2. Исправленный метод Эйлера. В исправленном методе Эйлера мы находим средний tg угла наклона касательной для двух точек: [I]x[m],y[m] и x[m+1],y[m]+h*y'[m][/I]. Соотношения, описывающие данный метод, имеют вид: y[m+1] = y[m] + h*F(x[m],y[m],h) (4) F(x[m],y[m],h) = ( y'[m] + f(x[m]+h, y[m]+h*y'[m]) )/2 (5) y'[m] = f(x[m],y[m]) (6) Исправленный метод Эйлера согласуется с разложением в ряд Тейлора вплоть до членов степени h^2, являясь, т.о. методом Рунге-Кутта второго порядка. 1.3. Модифицированный метод Эйлера. В данном методе мы находим tg угла наклона касательной в точке: x = x[m] + h/2; y = y[m] + (h/2)*y'[m] Соотношения, описывающие модифицированный метод Эйлера имеют вид: y[m+1] = y[m] + h*F(x[m],y[m],h) (7) F(x[m],y[m],h) = f( x[m]+h/2, y[m]+(h/2)*y'[m] ) (8) y'[m] = f(x[m],y[m]) (9) Модифицированный метод Эйлера также согласуется с разложением в ряд Тейлора вплоть до членов степени h^2, и также является методом Рунге-Кутта второго порядка. 1.4. Метод Рунге-Кутта четвертого порядка. Данный метод является одним из самых употребительных методов интегрирования дифференциальных уравнений. Этот метод применяется настолько широко, что в литературе его просто называют "методом Рунге-Кутта" без всяких указаний на его тип или порядок. Этот классический метод описывается системой следующих пяти соотношений: y[m+1] = y[m] + h*(k[1] + 2*k[2] + 2*k[3] + k[4])/6 (10) k[1] = f(x[m], y[m]) (11) k[2] = f(x[m]+h/2, y[m]+h*k[1]/2) (12) k[3] = f(x[m]+h/2, y[m]+h*k[2]/2) (13) k[4] = f(x[m]+h, y[m]+h*k[3]) (14) Ошибка приближения для этого метода равна e[t]=k*h^5. Заметим, что при использовании этого метода функцию необходимо вычислять четыре раза. Один из серьезных недостатков методов Рунге-Кутта состоит в отсутствии простых способов оценки их ошибки. 2. Методы прогноза и коррекции. Отличительной чертой методов Рунге-Кутта является то, что при вычислении следующей точки (x[m+1],y[m+1]) используется информация только об одной предыдущей точке (x[m],y[m]), но не нескольких. Кроме того, для методов Рунге-Кутта отсутствует достаточно простые способы оценки ошибки, что приводит к необходимости рассмотрения некоторых дополнительных методов решения ДУ. Отличительное свойство этих методов состоит в том, что с их помощью нельзя начать решение уравнения т.к. в них необходимо использовать информацию о предыдущих точках решения. Чтобы начать решение, имея только одну точку, определяемую начальными условиями, или для того, чтобы изменить шаг (h), необходим метод типа Рунге-Кутта. Поэтому приходится использовать разумное сочетание этих двух методов. Методы, которые мы рассмотрим, известны под общим названием методов прогно- за и корректировки. Как ясно из названия вначале "предсказывается" значение y[m+1], а затем используется тот или иной метод его "корректировки". Среди множества возможных формул прогноза и коррекции выберем по одному примеру, применимому ко многим практическим задачам. 2.1. Метод Адамса-Бошфора. Для прогноза используем формулу второго порядка y^(0)_[m+1] = y[m-1] + 2*h*f(x[m],y[m]) (15) где (0) - означает исходное приближение y[m+1], т.е. предсказанное значение. Непосредственно из написанной формулы следует, что с ее помощью нельзя вычис- лить y[1], т.к. для его вычисления потребовалась бы точка, расположенная перед начальной точкой y[0]. Чтобы начать решение с помошью метода прогноза и коррек- ции, для нахождения y[1] необходимо использовать метод типа Рунге-Кутта. Для коррекции возьмем формулу, похожую на исправленный метод Эйлера (4)-(6): y^(i)_[m+1] = y[m] + h(f(x[m],y[m])+f(x[m+1],y^(i-1)_[m+1]))/2 (16) для i=1,2,3, ... Итерационный процесс прекращается, когда | y^(i+1)_[m+1] - y^(i)_[m+1] | < Eps (17) для некоторого Eps>0. { Writeln('Labinskiy Nikolay aka e-moe (c) 2005'); Writeln('See you later on forum.sources.ru ;)'); } {$E+,N+,X+} const eps = 0.00001; x0 = 3; y0 = 3; type TFunc = function(x,y: extended): extended; function func(x,y: extended): extended; far; begin func := y*cos(x) - 2*sin(2*x); end; function Euler(f: TFunc; x0_,x_end,y0_: extended; n: word): extended; { где x0_ и y0_ - начальное условие x_end - точка, в которой необходимо вычислить результат n - количество шагов для вычисления результата } var i : word; { счетчик цикла } x,h : extended; { текущая точка и длина шага } res : extended; { переменная для накопления конечного результата функции} begin h:= (x_end - x0_)/n; { Находим длину шага } res:= y0_; { устанавливаем начальные значения} x:=x0_; for i:=1 to n do begin { вычисляем результат по методу Эйлера } res:=res+h*f(x,res); x:=x+h; { переходим к следующей точке } end; Euler:=res; { присваиваем конечный результат функции } end; function Euler2(f: TFunc; x0_,x_end,y0_: extended; n: word): extended; { где x0_ и y0_ - начальное условие x_end - точка, в которой необходимо вычислить результат n - количество шагов для вычисления результата } var i : word; { счетчик цикла } x,h : extended; { текущая точка и длина шага } res : extended; { переменная для накопления конечного результата функции} begin h:= (x_end - x0_)/n; { Находим длину шага } res:= y0_; { устанавливаем начальные значения} x:=x0_; for i:=1 to n do begin { вычисляем результат по исправленному методу Эйлера } res:=res+h*(f(x,res)+f(x+h,res+h*f(x,res)))/2; x:=x+h; { переходим к следующей точке } end; Euler2:=res; { присваиваем конечный результат функции } end; function Euler3(f: TFunc; x0_,x_end,y0_: extended; n: word): extended; { где x0_ и y0_ - начальное условие x_end - точка, в которой необходимо вычислить результат n - количество шагов для вычисления результата } var i : word; { счетчик цикла } x,h : extended; { текущая точка и длина шага } res : extended; { переменная для накопления конечного результата функции} begin h:= (x_end - x0_)/n; { Находим длину шага } res:= y0_; { устанавливаем начальные значения} x:=x0_; for i:=1 to n do begin { вычисляем результат по модифицированному методу Эйлера } res:=res+h*f(x+h/2,res+(h/2)*f(x,res)); x:=x+h; { переходим к следующей точке } end; Euler3:=res; { присваиваем конечный результат функции } end; function RungeKutt(f: TFunc; x0_,x_end,y0_: extended; n: word): extended; { где x0_ и y0_ - начальное условие x_end - точка, в которой необходимо вычислить результат n - количество шагов для вычисления результата } var i : word; { счетчик цикла } x,h : extended; { текущая точка и длина шага } res : extended; { переменная для накопления конечного результата функции } k1,k2,k3,k4: extended; { вспомогательные переменные вычисления результата } begin h:= (x_end - x0_)/n; { Находим длину шага } res:= y0_; { устанавливаем начальные значения} x:=x0_; for i:=1 to n do begin { вычисляем результат по методу Рунге-Кутта 4го порядка } k1:=f(x,res); k2:=f(x+h/2,res+h*k1/2); k3:=f(x+h/2,res+h*k2/2); k4:=f(x+h,res+h*k3); res:=res+h*(k1+2*k2+2*k3+k4)/6; x:=x+h; { переходим к следующей точке } end; RungeKutt:=res; { присваиваем конечный результат функции } end; function AdmsBoshf(f: TFunc; x0_,x_end,y0_: extended; n: word): extended; { где x0_ и y0_ - начальное условие x_end - точка, в которой необходимо вычислить результат n - количество шагов для вычисления результата } var i : word; { счетчик цикла } x,h : extended; { текущая точка и длина шага } y1,y2 : extended; { переменные для вычисления следующей точки } tmp : extended; { временная переменная для уточнения результата } begin h:= (x_end - x0_)/n; { Находим длину шага } x:=x0_; { устанавливаем текущую точку } y1:=y0_+h*f(x+h/2,y0_+(h/2)*f(x,y0_)); { находим начальное значение модифицированным методом Эйлера } repeat { уточняем(корректируем) начальное значение } tmp:=y1; y1:=y0_+h*(f(x,y0_)+f(x+h,y1))/2; until abs(y1-tmp) < eps; x:=x+h; { переходим к следующей точке } y2:=y0_+2*h*f(x,y1); { прогнозируенм следующее значение } repeat { уточнаяем наш прогноз } tmp:=y2; y2:=y1+h*(f(x,y1)+f(x+h,y2))/2; until abs(y2-tmp) < eps; for i:=3 to n do { начинаем счет с 3 т.к. две точки уже найдены } begin x0_:=x; { переходим к следующей точке } x:=x+h; y0_:=y1; y1:=y2; y2:=y0_+2*h*f(x,y1); { прогнозируенм следующее значение } repeat { уточнаяем наш прогноз } tmp:=y2; y2:=y1+h*(f(x,y1)+f(x+h,y2))/2; until abs(y2-tmp) < eps; end; AdmsBoshf:=y2; { присваиваем конечный результат функции } end; begin writeln('Численное решение дифференциальных уравнений:'); writeln(#10,' y'' = y*cos(x) - 2*sin(2*x); y(0)=3; x_end=',(5*x0+3.5):5:5); writeln(#10,'Метод Эйлера:'); writeln('n=5 000, result: ',Euler(func,x0,5*x0+3.5,y0,5000):5:5); writeln('n=10 000, result: ',Euler(func,x0,5*x0+3.5,y0,10000):5:5); writeln('n=25 000, result: ',Euler(func,x0,5*x0+3.5,y0,25000):5:5); writeln(#10,'Исправленный метод Эйлера:'); writeln('n=50, result: ',Euler2(func,x0,5*x0+3.5,y0,50):5:5); writeln('n=100, result: ',Euler2(func,x0,5*x0+3.5,y0,100):5:5); writeln('n=250, result: ',Euler2(func,x0,5*x0+3.5,y0,250):5:5); writeln(#10,'Модифицированный метод Эйлера:'); writeln('n=50, result: ',Euler3(func,x0,5*x0+3.5,y0,50):5:5); writeln('n=100, result: ',Euler3(func,x0,5*x0+3.5,y0,100):5:5); writeln('n=250, result: ',Euler3(func,x0,5*x0+3.5,y0,250):5:5); writeln(#10,'Метод Рунге-Кутта:'); writeln('n=50, result: ',RungeKutt(func,x0,5*x0+3.5,y0,50):5:5); writeln('n=100, result: ',RungeKutt(func,x0,5*x0+3.5,y0,100):5:5); writeln('n=250, result: ',RungeKutt(func,x0,5*x0+3.5,y0,250):5:5); write(#10,'Press Enter to continue.'); readln; writeln(#10,'Метод Адамса-Бошфора:'); writeln('n=50, result: ',AdmsBoshf(func,x0,5*x0+3.5,y0,50):5:5); writeln('n=100, result: ',AdmsBoshf(func,x0,5*x0+3.5,y0,100):5:5); writeln('n=250, result: ',AdmsBoshf(func,x0,5*x0+3.5,y0,250):5:5); end. Автор статьи: Labinskiy Nikolay aka e-moe © 2005 Это сообщение было перенесено сюда или объединено из темы "Численное решение ДУ" |
Сообщ.
#5
,
|
|
|
Вычисление приближённого значения функций с помощью ряда Тейлора (Разложение функций в степенные ряды) (применяется также для нахождения пределов вообще)
К примеру, sin = x -(x3/3!)+(x5/5!)- ... Здесь описано как вычислять Sin, Cos, ArcTg, Exp. Другие функции пишутся аналогично. Значение функции вычисляется в процедуре FunctionValue. Только подставьте нужное название функции. Для того, чтобы вычислить приближённое значение по формуле, мы должны взять бесконечно большое число, чтобы приблизиться к истинному значению функции. Т.е. при n->infinity (бесконечность), получаем истинное значение функции в данной точке. На практике мы этого не станем делать, так как в Паскале имеются типы с ограниченным размером, такие как: Real, Single, Double, Extended. Что же теперь делать? Введём понятие точность вычисления и назовём эту переменную eps (epsilon). Точность вычисления может быть разная: 0.05, 0.001 и намного большая, 10-8. В большинстве случаев хватает точности 3 знака после нуля, поэтому возьмём eps=0.001, или в краткой форме записи для Паскаля, 1e-3, что означает 1*10-3. Типа Single для этой точности вполне хватит, на нём и остановимся. Если истинное значение функции в заданной точке равно, допустим, 1.2328, а вычисленное значение равно 1.2323, то при нашей точности вычисления разница по модулю между ними составляет 0.0005<eps, то есть ответ нас устраивает. Для вычисления бесконечного ряда нужно применить цикл, в котором мы каждый раз будем проверять, насколько точно вычислено требуемое значение X, то есть проверим:|Xnew - Xold| < eps ? где Xnew - это новое вычисленное значение, Xold - предыдущее вычисленное значение Как только в ходе вычисления внутри цикла мы выясним, что достигнута требуемая точность (eps), выходим из него. Дальше вычислять уже нет необходимости. Далее. Как вычислять эти значения? Смотрим на общую формулу: sin (x) = Sum [(-1)n * x2n+1 / (2n+1)!] , где n->infinity. Из чего состоит данная формула? Из вычисления степени числа (функция Power), факториала числа (функция Factorial) и вычисления знака (функция Sign). Для повторения вычисления этих функций смотрите прикреплённые вверху ссылки. Теперь надо собрать все эти функции в одном выражении, которое будет вычисляться в отдельной функции. { вычисление значения элемента с индексом [I]i[/I] ряда Тейлора для синуса} function Sine (x: single; i: integer): single; var k: integer; begin K := 2 * i + 1; Sine := ( Sign (i) * Power (x, k) / factorial (k) ); end; Что осталось? Объединить всё написанное в одной функции вычисления приближённого значения FunctionValue. function FunctionValue (x: single): single; А теперь программа: const eps: single = 1e-3; { Epsilon - необходимая точность, 1*10^(-3) } function Factorial (N: word): single; { N! } var f: single; begin Factorial := 1.0; if n = 0 then exit; f := 1.0; for n := 1 to n do f := f * n; Factorial := f; end; function Sign (p: integer): single; { (-1)^p } begin if Not Odd(P) then Sign := 1 else Sign := -1; end; function Power (base: single; N: integer): single; { степень N по основанию base } var k: integer; P: single; begin P := 1.0; for k := 1 to N do P := P * base; Power := P; end; {--------------------------------------------------------------------} function Sine (x: single; i: integer): single; var k: integer; begin K := 2 * i + 1; Sine := ( Sign (i) * Power (x, k) / factorial (k) ); end; function CoSine (x: single; i: integer): single; var k: integer; begin K := 2 * i; CoSine := ( Sign (i) * Power (x, k) / factorial (k) ); end; function ArcTg (x: single; i: integer): single; var k: integer; begin K := 2 * i + 1; ArcTg := ( Sign (i) * Power (x, k) / k ); end; function Exp (x: single; i: integer): single; begin Exp := Power (x, i) / factorial (i); end; {--------------------------------------------------------------------} function FunctionValue (x: single): single; var sum,old: single; index: integer; begin if x <= 0.0 then begin FunctionValue := 0.0; exit; end; index := 0; sum := 0.0; repeat { повторять } old := sum; { old - прежнее вычисленное значение функции } sum := sum + ArcTg (x, index); { Sin, Cos, Exp } index := index + 1; until abs (sum - old) < eps; { пока не достигнута необходимая точность } FunctionValue := sum end; begin writeln ( 'ArcTg (x) =', FunctionValue ( 1.0 / sqrt (3.0) ) : 8 : 3 ); { ArcTg } { writeln ( 'Cosinus (x) =', FunctionValue ( pi / 6.0 ) : 8 : 3 ); { CoSine } readln; end. Прикреплённый файлtailor.pas.zip (0.89 Кбайт, скачиваний: 619) |