Применение пакета MATHCAD к решению вычислительных задач
.pdf
1.5. Решение СЛАУ
Стандартным средством решения линейных систем является встроенная функция lsolve(A,b), которая решает систему на основе LU-разложения матрицы (с частичным выбором главного элемента). Ее аргументами являются квадратная невырожденная матрица A и вектор правых частей b.
2x1 + x2 −3x3 |
|
= 4 |
|||
Пример. Решим систему |
x1 + 2x2 + x3 |
|
= 5 |
||
|
3x −2x |
2 |
+ 2x |
3 |
= −1 |
|
1 |
|
|
||
|
2 |
1 |
−3 |
||
Зададим матрицу системы: |
A := |
1 |
2 |
1 |
|
|
|
|
|
|
|
|
|
3 |
−2 |
2 |
|
4
Зададим вектор правых частей: b := 5
−1
|
|
1 |
|
|
Найдем решение: |
x := lsolve(A,b) |
x = |
2 |
|
|
|
|
|
|
|
|
|
0 |
|
4
Для проверки умножим матрицу A на найденное решение A x = 5
−1
Произведение совпало с правой частью b.
1.6.Интерполяция функций
Впакете Mathcad имеется возможность интерполяции кубическими сплайнами. Для этого предназначена пара функций cspline(…) – interp(…). Если заданы массивы узлов интерполяции X и соответствующих значений функции Y, то функция cspline(X,Y) возвращает массив коэффициентов для построения кубического сплайна. Этот массив должен быть первым аргументом функции interp(S,X,Y,t). Два следующих аргумента – это массивы исходных данных,
11
последний аргумент – точка, в которой вычисляется значение построенного интерполяционного сплайна.
Пример. Построим интерполяционный сплайн для функции, заданной
таблицей значений |
x |
|
−1 |
|
−0.5 |
|
0 |
|
0.5 |
|
1 |
|
|
|
|
|
|
||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
y |
|
3 |
|
2 |
|
1 |
|
2.5 |
|
4 |
||
|
|
|
|
|
|
|||||||
Задаем исходные данные:
|
−1 |
|
|
3 |
|
|
−0.5 |
|
|
2 |
|
|
|
|
|
||
x := |
0 |
|
y := |
1 |
|
|
0.5 |
|
|
2.5 |
|
|
1 |
|
|
4 |
|
|
|
|
|
Получаем коэффициенты многочленов:
s := cspline(x,y)
Задаем интерполяционный сплайн:
g(t) := interp(s ,x,y ,t)
Строим график полученного сплайна и точечный график исходных данных:
g(t) |
4 |
|
|
y |
2 |
|
|
|
|
|
|
|
1 |
0 |
1 |
|
|
t,x |
|
Описанным образом можно построить кубический сплайн дефекта один с условиями отсутствия узла. Помимо этого, имеется возможность построения указанного сплайна с другими дополнительными условиями. А именно, использование функции pspline() для получения массива коэффициентов позволяет построить сплайн, у которого в граничных точках равна нулю третья производная. А функция lspline() дает указанный сплайн с равными нулю в граничных точках второй и третьей производными.
12
1.7. Решение задачи Коши
Одной из простейших встроенных функций для решения задачи Коши для
ОДУ 1 порядка вида y'= f (t, y) |
является функция rkfixed(y0,t0,tN,N,F), в |
y(t0) = y0 |
|
которой реализован метод Рунге-Кутты 4-го порядка. Ее первый аргумент y0 – начальное значение решения, второй и третий аргументы t0, tN – начальное и конечное значения аргумента (начало и конец отрезка, на котором решается задача), четвертый аргумент N – количество точек на указанном отрезке, в которых будет найдено решение и последний аргумент F – имя функции – правой части дифференциального уравнения. Функция rkfixed(y0,t0,tN,N,F) возвращает матрицу из двух столбцов и N+1 строки. В первом (по счету) столбце содержатся узлы, на которые разбит отрезок [t0, tN], во втором столбце – значения искомого решения в соответствующих узлах отрезка.
y'= sin(t y)
Пример. Решим задачу Коши y(0) =1
зададим начальные данные: |
t0 := 0 y0 := 1 |
зададим правую часть ОДУ: |
f(t ,y) := sin(t y) |
решим задачу на отрезке длины 2, разбитом на 10 шагов:
y := rkfixed(y0 ,t0 ,t0 + 2 ,10 ,f)
|
|
0 |
|
1 |
|
0 |
|
0 |
1 |
|
1 |
|
0.2 |
1.02 |
|
2 |
|
0.4 |
1.082 |
|
3 |
|
0.6 |
1.189 |
y = |
4 |
|
0.8 |
1.343 |
5 |
|
1 |
1.534 |
|
|
|
|||
|
6 |
|
1.2 |
1.727 |
|
7 |
|
1.4 |
1.868 |
|
8 |
|
1.6 |
1.923 |
|
9 |
|
1.8 |
1.9 |
|
10 |
|
2 |
1.822 |
|
|
|
|
|
Выделим из полученного результата массив аргументов (вырежем из матрицы y нулевой (первый по счету)
X := y 0
13
Выделим из полученного результата массив значений решения (вырежем из матрицы y первый (второй по счету)
Y := y 1
Таким образом, получено решение:
|
|
0 |
|
|
|
|
|
|
|
|
0 |
||
|
0 |
0 |
|
|
|
|
|
|
0 |
1 |
|||
|
1 |
0.2 |
|
|
|
|
|
|
1 |
1.02 |
|||
|
2 |
0.4 |
|
2 |
1.082 |
|
|
3 |
0.6 |
|
3 |
1.189 |
|
X = |
4 |
0.8 |
Y = |
4 |
1.343 |
|
5 |
1 |
5 |
1.534 |
|||
|
|
|||||
|
6 |
1.2 |
|
6 |
1.727 |
|
|
7 |
1.4 |
|
7 |
1.868 |
|
|
8 |
1.6 |
|
8 |
1.923 |
|
|
9 |
1.8 |
|
9 |
1.9 |
|
|
|
|
|
10 |
1.822 |
|
|
10 |
2 |
||||
|
|
|
|
|
|
построим график найденного решения:
|
2 |
|
|
|
Y |
1.5 |
|
|
|
|
1 |
0 |
1 |
2 |
|
|
|
X |
|
Определим, например, значение y(1).
Для этого найдем в массиве X соответствующее значение аргумента. Это X5 = 1 Теперь возьмем из массива Y
элемент с тем же номером: Y5 = 1.534 Это и есть искомое значение.
14
Ту же самую функцию rkfixed(y0,t0,tN,N,F) можно использовать при решении задачи Коши для систем ОДУ 1-го порядка. В этом случае ее первый аргумент y0 должен быть вектором начальных приближений, а последний аргумент F(t,y) должен быть вектор-функцией правых частей (зависящей от скаляра t и вектора y). Возвращать же функция rkfixed(y0,t0,tN,N,F) будет матрицу с большим количеством столбцов. В первом (по счету) столбце опять будут содержаться узлы сетки, а в остальных – значения искомых функций в этих узлах (и количество таких столбцов будет совпадать с порядком системы).
Пример. Решим задачу Коши для системы
u' = sin(t v)
v′ = u
u(0) =1, v(0) = 0.5
зададим начальные данные: |
t0 := 0 |
y0 := |
|
1 |
|
|
|
0.5 |
|
||||
|
|
|
||||
зададим правую часть ОДУ: |
sin(t y1) |
|
||||
f(t ,y) := |
y0 |
|
|
|
||
|
|
|
|
|
||
решим задачу на отрезке длины 2, разбитом на 10 шагов:
y := rkfixed(y0 ,t0 ,t0 + 2 ,10 ,f)
|
|
0 |
|
1 |
2 |
|
0 |
|
0 |
1 |
0.5 |
|
1 |
|
0.2 |
1.02 |
0.611 |
|
2 |
|
0.4 |
1.082 |
0.746 |
|
3 |
|
0.6 |
1.189 |
0.911 |
y = |
4 |
|
0.8 |
1.343 |
1.113 |
5 |
|
1 |
1.534 |
1.359 |
|
|
|
||||
|
6 |
|
1.2 |
1.727 |
1.66 |
|
7 |
|
1.4 |
1.868 |
2.028 |
|
8 |
|
1.6 |
1.923 |
2.476 |
|
9 |
|
1.8 |
1.9 |
3.025 |
|
10 |
|
2 |
1.822 |
3.694 |
|
|
|
|
|
|
Выделим из полученного результата массив аргументов: X := y 0
Выделим из полученного результата первую функцию: U := y 1
15
Выделим из полученного результата вторую функцию: V := y 2
Построим графики найденного решения: |
|
Первая функция: |
Вторая функция: |
|
4 |
|
|
|
4 |
|
|
U |
2 |
|
|
V |
2 |
|
|
|
0 0 |
1 |
2 |
|
0 0 |
1 |
2 |
|
|
X |
|
|
|
X |
|
1.8. Блок Given -- Find
Достаточно универсальным средством решения различных вычислительных задач является блок Given – Find. Покажем, как пользоваться данным средством на примере решения систем нелинейных уравнений. (К тому же, других способов решения таких систем Mathcad не имеет.) Между ключевыми словами записываются все уравнения системы, которую необходимо решить. При этом для знака равенства следует использовать «логическое равно» (набирается с помощью клавиш [Ctrl]= или выбирается с панели логических символов). После ключевого слова find в скобках указываются искомые величины.
Пример. Найдем точки пересечения эллипса x2 + 2y2 |
|
6 и |
|
||
|
прямой x − 2y 
0
задаем исходные данные: given
x2 + 2y2 |
|
6 |
||
|
||||
|
||||
x − 2y |
|
|
0 |
|
|
|
|
||
|
|
|||
вызываем функцию решения системы:
16
find(x,y) → |
2 |
−2 |
|
1 |
−1 |
||
|
В итоге получаем, что данная задача имеет два решения (они расположены в выданном ответе по столбцам):
x1 = 2 y1 = 1 и x2 = −2 y2 = −1
Для проверки построим графики заданных кривых. Выражая y из каждого уравнения, получаем:
для эллипса y0(x) := |
6 − x2 |
, |
y1(x) := − |
6 − x2 |
|
|
|
|
2 |
|
|
2 |
|
|
|
для прямой y2(x) := |
x |
|
|
|
|
|
|
Теперь строим график |
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
y0(x) |
|
|
|
1 |
|
|
|
y1(x) |
|
|
|
|
|
|
|
y2(x) |
3 |
2 |
1 |
0 |
1 |
2 |
3 |
|
|
|
1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
x |
|
|
|
По графику хорошо видно, что решений действительно два и они найдены правильно.
17
Часть II
Встроенный язык программирования
Пакет Mathcad обладает некоторым арсеналом средств, позволяющим создавать пользовательские функции и реализовывать таким образом практически любые алгоритмы вычислительной математики. Все соответствующие средства собраны на панели программирования (programming toolbar). Следует сразу подчеркнуть, что в текст программы все операторы, ключевые слова и т.д. должны вставляться только с панели программирования и ни в коем случае не вводиться посимвольно с клавиатуры. Никаких принципиальных отличий в использовании операторов во встроенном языке пакета Mathcad нет. Ниже для удобства восприятия все алгоритмы продублированы на языке Паскаль (предполагается, что он уже освоен читателем). С одной стороны, это позволяет лучше понять структуру и детали алгоритма, а с другой – увидеть те небольшие отличия, которые есть между двумя языками программирования.
Любая программа состоит из заголовка и тела программы. Заголовок представляет собой название программы (произвольная последовательность стандартно допустимых символов) и параметров. Параметры программы задаются простым перечислением и не требуют явного указания их типа. Тип распознается автоматически по использованию того или иного параметра в теле программы. Это могут быть числа (целые, вещественные, комплексные), вектора (массивы), матрицы, функции и т.д. За заголовком следует оператор присваивания и начинается собственно тело (текст) программы. Оно ограничено жирной вертикальной чертой, которая создается по нажатию на кнопку Add Line и служит аналогом операторных скобок (begin ... end в Паскале, { ... } в Си и т.д.) Если необходимо, чтобы программа после завершения своей работы возвращала какое-либо значение, то имя переменной, в которой это значение хранится, необходимо указать в самой последней строке программы. Внутренние переменные в программе также не требуют предварительного описания с указанием типа данных. Фактически, определением переменной является первый оператор присваивания, в левой части которого она встречается.
Теперь разберем на примерах основные конструкции встроенного языка.
Начнем с оператора присваивания ( 
← 
) и цикла for.
Для начала напишем простую программу, вычисляющую норму вектора,
n
т.е. величину 
x
1 = ∑xi (по сути дела, найдем сумму модулей n чисел). При
i=0
этом не забываем, что нумерация массивов в пакете Mathcad по умолчанию ведется с нуля.
18
Программа 1. |
Вычисление |
|
|
|
x |
|
|
|
1 |
|
|
|
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||||||
type vector = array[0..100] of double; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||
function norma1(x: vector; n: word): |
|
|
|
|
|
|
normа1(x,n) := |
|
s ← |
|
x0 |
|
|
|
|
|
||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||||||||
double; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
var i: word; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
for |
i 1.. n |
||||||||
s: double; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
s ← s + |
|
xi |
|
||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||
begin |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
s := x[0]; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
s |
|
|
|
|
|
|
|
||
for i := 1 to n do |
s := s + abs(x[i]); |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||
normа1 := s; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
end; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||
|
Пример вызова программы 1. |
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|||
|
Зададим вектор x: |
|
x := |
−3 |
|
|
|
|
|
|
|
|
||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1 |
|
|
|
|
|
|
|
|
||
|
Так как нумерация начинается с нуля, то |
|
|
|
|
|
|
|
||||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
n := 2 |
|
|
|
|
|
|
|
|
|
|
|
|||||||||
|
(номер последнего элемента равен двум) |
|
|
|
|
|
|
|
||||||||||||||||||||||||
|
Вызовем программу: |
|
normа1(x,n) = 6 |
|
|
|
|
|
|
|
||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||||||||
Ниже приведена еще одна простая программа, вычисляющая другую |
||||||||||||||||||||||||||||||||
норму вектора |
|
x |
|
|
|
∞ |
= max |
|
xi |
|
(теперь ищется максимальное по модулю число |
|||||||||||||||||||||
|
|
|
|
|
||||||||||||||||||||||||||||
|
|
|
|
|
|
|
|
0≤i≤n |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
среди заданных).
На этом примере хорошо видна особенность условного оператора во встроенном языке пакета Mathcad: перед ключевым словом if пишется выполняемый оператор (или последовательность операторов), а после ключевого слова – условие, в случае истинности которого эти операторы будут выполнены.
19
|
Программа 2. |
Вычисление |
|
x |
|
|
|
∞ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
type |
vector = array[0..100] of double; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
function normai(x: vector; n: word): |
|
normai(x,n) := |
|
|
|
s ← |
|
x0 |
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
|
|
|
|||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||
double; |
|
|
|
|
|
|
|
|
|
|
|
|
|
for |
i 1.. n |
|
||||||||||||
var |
i: word; s: double; |
|
|
|
|
|
|
|
|
|
|
|
|
s ← |
|
xi |
|
if s < |
|
xi |
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||
begin |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||||
s := x[0]; |
|
|
|
|
|
|
|
|
|
|
|
|
|
s |
|
|
|
|
|
|
|
|
|
|
|
|
||
for i := 1 to n do |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
if s < abs(x[i]) then |
s := abs(x[i]); |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
normаi := s; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
end; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
Пример вызова программы 2. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|||||
|
|
Зададим тот же вектор x: |
|
|
|
|
x := |
−3 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
|
1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
Так как нумерация начинается с нуля, то |
|
|
|
|
n := 2 |
|
|
|||||||||||||||||||
|
|
Вызовем программу: |
|
normai(x,n) = 3 |
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Теперь немного усложним задачу: будем искать максимум не среди элементов массива, а на множестве значений функции f(x) на заданном отрезке [a,b]. Для этого зададим на отрезке [a,b] множество точек, расположенных равномерно с шагом h, и найдем максимальное значение функции в этих точках. Необходимо, конечно, понимать, что такой алгоритм не ищет точное значение маскимума функции. Он возвращает лишь некоторое приближение к нему, точность которого тем выше, чем большее количество точек мы рассмотрим на отрезке [a,b]. В приведенной ниже программе количество точек фиксировано (и равно 100). Для увеличения точности следует увеличить этот параметр в тексте программы (при вычислении величины h), а еще лучше – ввести дополнительный аргумент, задающий количество рассматриваемых на отрезке
[a,b] точек.
20