Добавил:

Hidi Опубликованный материал нарушает ваши авторские права? Сообщите нам.

Вуз:

Национальный исследовательский университет «МЭИ»

Предмет:

MathCad/MatLab/Maple

Файл:

Применение пакета MATHCAD к решению вычислительных задач

.pdf

Скачиваний:

253

Добавлен:

26.05.2014

Размер:

379.6 Кб

Скачать

☆

<<< < Предыдущая 1 23 / 33

Программа 3.

Поиск максимума функции на отрезке.

type

func = function(t: double):

double;

function maxf(f: func; a,b: double):

maxf(f ,a ,b) :=

x ← a

double;

s ← f(a)

var

s,x,h: double;

begin

h ← b − a

x := a;

100

s := f(a);

while x < b

h := (b - a) / 100;

while x < b do

x ← x + h

begin

s ← f(x) if s < f(x)

x := x + h;

if s < f(x) then s := f(x);

end;

maxf := s;

end;

Пример вызова программы 3.

Найдем максимум модуля

производной функции g(x) =

sin 2 (x)

на отрезке [1, 1.5]

Определим нужную функцию:

sin(x)2

g(x) :=

+ 1

Определим модуль ее производной:

h(x) :=

g(x)

Зададим границы отрезка:

a := 1

b := 1.5

Найдем максимум:

maxf(h ,a ,b) = 0.227

Для проверки построим график:

0.2

h(x)

0.1

1.1

1.2

1.3

1.4

1.5

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

h(1) = 0.227

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

Следующая программа уже может быть использована для решения серьезной задачи – поиска корня уравнения вида x =ϕ(x) методом простых итераций. При этом предполагается, что необходимая для работы метода величина q = maxϕ'(x) уже найдена. Напомним, что метод простой итерации

заключается в последовательном вычислении приближений xn+1 =ϕ(xn ) до тех пор, пока не выполнится условие xn+1 − xn ≤1−qq ε , гарантирующее нахождение

решения с точностью ε .

В качестве ответа приведенная ниже программа выдает массив из двух значений: корня уравнения и количества итераций. Входными параметрами являются функция - правая часть уравнения ϕ(x) (она обозначена через g), знаменатель метода q, начальное приближение x0 и точность e. Две переменные x и y в теле цикла отвечают за предыдущее и последующее приближения к корню соответственно, а их разность позволяет следить за выполнением критерия окончания.

Программа 4.

Метод простой итерации.

type

func = function(t: double):

double;

vect2 = array[0..1] of double;

siter(g,q ,x0

,e) :=

x ← x0

function siter(f: func; q, x0, e: double):

y ← g(x)

vect2;

k ← 1

var

x,y: double;

1 − q

k: word;

e1 ←

begin

x := x0;

while

x − y

> e1

y := g(x);

x ← y

k := 1;

e1 := ((1 - q) / q ) * e;

y ← g(x)

while abs(x - y) > e1 do

k ← k + 1

begin

x := y;

res ←

y := g(x);

k := k + 1;

res

end;

res[0] := y;

res[1] := k;

siter := res;

end;

Пример вызова программы 4. Найдем корень уравнения				x −	sin2	(x)	−1	= 0
Пример вызова программы 4. Найдем корень уравнения				x −	4		−1	= 0
					4
Зададим функцию:		2
f(x) := x −	sin(x)	2	− 1

	4
	4

Оставляя x слева и перенося остальные слагаемые вправо, получим уравнение вида x = g(x) . Зададим полученную функцию g(x):

g(x) := (sin(x))2 + 1 4

Построим график заданной функции f(x):

0.5

f(x)

1.2

1.4

1.6

1.8

0.5

По графику видно, что в качестве начального приближения можно выбрать значение 1.2, так как примерно там происходит пересечение графиком оси OX.

Определим знаменатель метода q = max g'(x) . (Это подробно расписано в

предыдущем примере). По графику функции f(x) видно, корень расположен на отрезке [1, 1.4] и, следовательно, можно ограничиться рассмотрением только этого отрезка.

h(x) :=		d	g(x)	q := maxf(h ,1 ,1.4)
		dx

Получили значение	q = 0.227

Оно меньше единицы, следовательно метод должен сходиться.

Зададим начальное приближение и точность:

x0 := 1.2

e := 10−8

Вызовем функцию, реализующую метод простой итерации:

y := siter(g,q ,x0,e)		1.221
Получим следующий результат:	y =
		9

Если нам нужно только значение корня, то вычленим его:			X := y0
Таким образом, ответом (корнем) является величина			X = 1.22056858

(Поскольку корень был найден с точностью 10-8, то необходимо настроить формат вывода числа Х так, чтобы отображались все 8 верных знаков после десятичной точки.)

Внесем в написанную программу небольшое изменение: введем предельное количество итераций, по достижении которого программа будет прерывать свою работу. (Это может пригодиться, например, в том случае, когда знаменатель q определен неправильно, и метод может расходиться.) Используем для этих целей оператор break, который будет срабатывать при достижении счетчиком числа итераций значения 10000. В такой ситуации произойдет выход из цикла while и начнут выполняться следующие за ним операторы. Кроме того, программа должна возвращать информацию о том, по какой причине закончилась ее работа: по выполнению критерия окончания или по достижению предельного числа итераций. Для этого введем дополнительную переменную key, которая будет равна 1 при нормальной работе метода и становиться нулем при превышении предельного числа итераций. Соответственно, в массив возвращаемых значений необходимо добавить эту переменную. Во всем остальном новая программа работает так же, как предыдущая.

Программа 5. Метод простой итерации с прерыванием.

type

func = function(t: double): double;

siter(g,q ,x0,e) :=

x ← x0

vect2 = array[0..1] of double;

y ← g(x)

function siter(f: func; q, x0, e: double):

k ← 1

vect2;

key ← 1

var

x,y: double;

1 − q

k: word;

e1 ←

begin

x := x0;

y := g(x);

while

x − y

> e1

k := 1;

key := 1;

x ← y

e1 := ((1 - q) / q ) * e;

y ← g(x)

while abs(x - y) > e1 do

begin

k ← k + 1

x := y;

k > 10000

y := g(x);

k := k + 1;

break

if k > 10000

key ← 0

then begin key := 0; exit; end;

end;

siter[0] := y;

res ←

siter[1] := k;

siter[2] := key;

key

end;

res

Чтобы набрать новый вариант программы, необходимо после того, как прописан оператор if, нажать кнопку Add Line. Тогда под оператором if появится вертикальная линия с несколькими полями ввода.

Вызов программы 5 производится точно так же, как и программы 4, поэтому подробно он не прописан.

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

Программа 6.	Метод бисекции.
type		bisec(f ,a ,b ,e) :=	an ← a
func = function(t: double):		bisec(f ,a ,b ,e) :=	an ← a
double;			bn ← b
vect2 = array[0..1] of double;			k ← 0
			k ← 0
function bisec(f: func; a, b, e:			while (bn − an) > 2 e
double): vect2;				xn ←		an + bn
double): vect2;
var xn, an, bn: double;
var xn, an, bn: double;
k: word;				2
k: word;				bn ← xn			if f(an) f(xn) ≤ 0
begin				bn ← xn			if f(an) f(xn) ≤ 0
an := a; bn := b; k := 0;				an ← xn			otherwise
while bn - an > 2*e do				an ← xn			otherwise
while bn - an > 2*e do				k ← k + 1
begin
begin
xn := (an + bn) / 2;			xn ←		an + bn
if f(an)*f(xn) <= 0					an + bn
if f(an)*f(xn) <= 0
then bn := xn			2
then bn := xn					xn
else an := xn;					xn
k := k + 1;			res ← k
end;			res
xn := (an + bn) / 2;			res
xn := (an + bn) / 2;
bisec [0] := xn;
bisec [1] := k;
end;

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

точность, а возвращается массив из двух значений – корня функции и количества произведенных итераций. Внутренние переменные an и bn содержат границы текущего отрезка локализации, xn – середину этого отрезка. В том случае, когда знак функции на левом конце отрезка не совпадает с ее знаком в середине

отрезка ( f(an) f(xn) ≤ 0), в качестве нового отрезка локализации выбирается левая половина (т.е. правой границей отрезка становится предыдущая середина:

bn ← xn). В противном случае выбирается правая половина (т.е. левой границей отрезка становится предыдущая середина: an ← xn).

Пример вызова программы 6. Найдем корень уравнения из предыдущего примера.

Зададим функцию:	f(x) := x −	sin(x)	2	− 1

		4

Зададим точность:	e := 10−8

Вызовем метод бисекции на отрезке локализации [1, 1.4] (см. график в предыдущем примере):

y := bisec(f ,1 ,1.4 ,e)
Получим результат:	y =	1.22056858
			25

Напоследок, остается еще раз сказать о возвращаемых значениях. Во всех приведенных выше программах возвращается либо одно значение, либо несколько однотипных. Однако, иногда бывает необходимо возвращать данные разных типов. Технически это делается точно так же: заполняется массив значений и возвращается этот массив. Единственным, на что здесь следует обратить внимание, является работа с массивами. А именно, если одним из элементов возвращаемого массива данных является также массив, то в выдаваемом векторе будет отображаться не содержимое этого массива, а размерность (взятая в фигурные скобки). Чтобы увидеть сами элементы, необходимо произвести разыменование. Проиллюстрируем сказанное на примере программы, определяющей оптимальное значение параметра для метода простых итераций (для СЛАУ). Напомним, что оптимальным параметром

является величина a = m +2M , где m и M – минимальное и максимальное

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

Программа 7. Определение оптимального параметра для метода простых итераций.

В языке Pascal создание типов данных, аналогичных массиву res из программы Mathcad’а невозможно.

param(A) := z ← eigenvals(A) m ← min(z)

M ← max(z)

a ← m + M if m > 0

a ← 0 otherwise res ← za

res

Пример вызова программы 7.		Определим оптимальный параметр метода
				17		−5	6	−6
					−5	11	−6	4
простых итераций для матрицы		A =							.
					6	−6	9	−8

						4	−8	13
				−6
		17			−5	6	−6
			−5		11	−6	4
Зададим матрицу:
	A :=			6	−6	9	−8

					4	−8	13
		−6
Вызовем написанную выше программу:							y := param(A)

Получим результат:	y =	{4,1}
		0.062

Как видим, значение параметра выдается сразу: y1 = 0.062

А вместо собственных чисел выдается шифрованная запись {4,1}.

Это означает, что первым (по счету) элементом массива y является массив, состоящий из 4-х строк и 1-го столбца. Чтобы его увидеть, обратимся к соответствующему элементу массива y:

30.4

1.928 y0 = 9.825

7.847

Теперь все собственные числа видны.

Чтобы обратиться к какому-либо из них, необходимо ставить скобки:

(y0)2 = 9.825

Или вводить новую переменную:

z := y0

z2 = 9.825

Библиографический список

1.Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. М.: Изд-во МЭИ, 2003.

2.Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1999.

3.Амосова О.А., Григорьев В.П., Зайцева С.Б. Вычислительные методы с применением математического пакета MATHCAD. М.: Изд-во МЭИ, 2000.

4.Амосова О.А., Акимова И.Г., Зайцева С.Б. Сборник заданий к лабораторным занятиям по вычислительной математике. М.: Изд-во МЭИ, 2002.

5.Дьяконов В. Mathcad 8 / 2000: специальный справочник. СПб. Изд-во Питер, 2000.

<<< < Предыдущая 1 23 / 33

Соседние файлы в предмете MathCad/MatLab/Maple

#
26.05.2014126.92 Кб50Перевод единиц измерения на русский язык.xmcd
#
26.05.2014379.6 Кб253Применение пакета MATHCAD к решению вычислительных задач.pdf
#
26.05.201436.42 Кб52Системы Единиц.mcd