Лучший способ сделать итерационную схему - PullRequest
Новая
       6

Лучший способ сделать итерационную схему

3 голосов
/ 13 января 2012

Надеюсь, об этом раньше не спрашивали, если да, то приношу свои извинения.

РЕДАКТИРОВАТЬ: Для ясности будут использоваться следующие обозначения: прописные буквы жирным шрифтом для матриц, строчные буквы жирным шрифтом для векторов икурсив для скаляров.

Предположим, x0 - вектор, A и B - матричные функции, а f - это векторная функция.

Я ищу лучший способ сделать следующую итерационную схему в Mathematica: 

A0 = A(x0), B0=B(x0), f0 = f(x0)
x1 = Inverse(A0)(B0.x0 + f0)

A1 = A(x1), B1=B(x1), f1 = f(x1)
x2 = Inverse(A1)(B1.x1 + f1)

...

Я знаю, что for-loop может добиться цели, но я не совсем знаком с Mathematica, и я обеспокоен тем, что это самый эффективный способ сделать это.Это оправданная проблема, поскольку я хотел бы определить функцию u(N):=xN и использовать ее в дальнейших вычислениях.

Наверное, мои вопросы:

Какой самый эффективный способ программирования схемы?

Является ли RecurrenceTable подходом?

РЕДАКТИРОВАТЬ

Это было немного сложнее, чем я думал.Я предоставляю более подробную информацию, чтобы получить более подробный ответ.

Перед повторением у меня возникают проблемы с пониманием того, как программировать функции A , B и f .

Матрицы A и B являются функциями временного шага dt = 1 / T и шаг пробела dx = 1 / M , где T и M - количество точек в { 0 , 0 } регион.Это также верно для векторной функции f .

Зависимость A , B и f on x довольно сложно:

A и B являются верхней и нижней треугольной матрицами (как матрица трехдиагональная ; я полагаю, что мы можем назвать их мультидиагональные ) с определенными постоянными значениями на своих диагоналях.

При заданной точке 0 , мне нужно определить его репрезентативную xn в сетке (ближайшей), и затем заменить строку nth A и B с функцией v ( x ) (конечно, транспонированной) и строкой nth f с функцией w ( x ).

Суммируя, A = A ( dt , dx , xs , x ).То же самое верно для B и f .

Затем мне нужно сделать цикл, упомянутый выше, чтобы определить u ( x) = step[T].

Надеюсь, я все объяснил.

Ответы [ 2 ]

5 голосов
/ 13 января 2012

Я не уверен, что это лучший метод, но я бы просто использовал простое старое запоминание. Вы можете представить отдельный шаг как 

xstep[x_] := Inverse[A[x]](B[x].x + f[x])

, а затем 

u[0] = x0
u[n_] := u[n] = xstep[u[n-1]]

Если вы заранее знаете, сколько значений вам нужно, и по какой-то причине выгодно предварительно вычислить их все (например, вы хотите открыть файл, использовать его содержимое для вычисления xN, а затем освободить память), вы может использовать NestList. Вместо двух предыдущих строк вы должны сделать 

xlist = NestList[xstep, x0, 10];
u[n_] := xlist[[n]]

Это сломается, если n > 10, конечно (очевидно, измените 10 в соответствии с вашими фактическими требованиями).

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

2 голосов
/ 13 января 2012

Я, вероятно, написал бы функцию, которая принимает A0, B0, x0 и f0, а затем возвращает A1, B1, x1 и f1 - скажем, 

step[A0_?MatrixQ, B0_?MatrixQ, x0_?VectorQ, f0_?VectorQ] := Module[...]

Я бы тогда Nest эту функцию,Трудно быть более точным без более точной информации.

Кроме того, если ваша процедура числовая, вы, конечно, не хотите вычислять Inverse[A0], поскольку это не числовая стабильная операция.Вместо этого вы должны написать 

A0.x1 == B0.x0+f0

, а затем использовать численно устойчивый решатель, чтобы найти x1.Конечно, Mathematica LinearSolve предоставляет такой алгоритм.

...