Функция Python для решения Ax = b путем обратной подстановки - PullRequest
Новая
       26

Функция Python для решения Ax = b путем обратной подстановки

4 голосов
/ 22 сентября 2010

Хорошо, для моего класса численных методов у меня есть следующий вопрос:

Напишите функцию Python для решения Ax = b путем обратной замены, где A - верхняя треугольная невырожденная матрица.Код MATLAB для этого приведен на стр. 190, который вы можете использовать в качестве руководства для псевдокода, если хотите.Функция должна принимать в качестве входных данных A и b и возвращать x.Ваша функция не должна проверять, что A неособо.То есть, предположим, что в вашу функцию будет передан только неособой А.

Код MATLAB, на который он ссылается:фрагмент с верхней треугольной тестовой матрицей (не уверен, что она неособая! Как мне проверить единственность?): 

from scipy import mat
c=[3,2,1]
U=([[6,5,1],[0,1,7],[0,0,2]])
a=0
x=[]
while a<3:
    x.append(1)
    a=a+1

n=3
i=n-1
x[n-1]=c[n-1]/U[n-1][n-1]
while i>1: 
    x[i]=c[i]
    j=i+1
    while j<n-1:
        x[i]=x[i]-U[i][j]*x[j];
    x[i]=x[i]/U[i][i]
    i=i-1
print mat(x)

Ответ, который я получаю, - [[1 1 0]] для xЯ не уверен, правильно ли я это делаю.Я предполагаю, что это неправильно и не могу понять, что делать дальше.Любые подсказки?

Ответы [ 3 ]

2 голосов
/ 22 сентября 2010

Одна проблема, которую я вижу, состоит в том, что ваш ввод состоит из целых чисел, что означает, что Python собирается сделать целочисленное деление на них, что превратит 3/4 в 0, когда то, что вы хотите - это деление с плавающей запятой.Вы можете указать Python делать деление с плавающей запятой по умолчанию, добавив 

from __future__ import division

в начало вашего кода.Исходя из использования scipy, я предполагаю, что вы используете Python 2.x здесь.

2 голосов
/ 22 сентября 2010 
j=i+1
while j<n-1:
    x[i]=x[i]-U[i][j]*x[j];

бесконечно ... и никогда не исполняется

Ваша индексация fubared: 

for i in range(n-2,-1,-1):
....
    for j in range(i+1,n):

Обратите внимание, диапазон наполовину открыт в отличие от Matlab

1 голос
/ 22 сентября 2010

Вы спрашиваете, как проверить особенность верхней треугольной матрицы?

Пожалуйста, не вычисляйте детерминант!

Просто посмотрите на диагональные элементы.Какие из них равны нулю?Есть ли ноль?

Как насчет эффективной числовой сингулярности?Сравните наименьшее абсолютное значение с наибольшим абсолютным значением.Если это отношение меньше, чем что-то порядка eps, оно фактически является единичным.

...