Python-сообщество

shevzoom · Ноя. 5, 2020 22:44:35

Переписывал программу с R на Python, но почему-то значения y(x) на пайтоне слишком скачут в зависимости от N, при n=100 порядок ординат десятки, при 1000 порядок сотни, на R как было в пределах +- от 0 до 1, так и остается, получается какая - та неустойчивость решения на python. Кто-то может объяснить с чем мб связано дело?
Вот на R:

 require(ggplot2)
ua = 0.8
ub = 1
A = function(x){
  return( log(3 + sqrt(1 + x) ))
}
B = function(x){
  return( exp(-x / 2) )
}
C = function(x){
  return ( sqrt((3 - x) / 2) )
}
F = function(x){
  return ( 2 - abs(1 - x) )
}
run = function(N){
  h = 2.4/(N + 1)
  
  an = function(x){
    return( -(2 * A(x) + B(x) * h) )
  }
  bn = function(x){
    return( 4 * A(x) + 2 * C(x) *h^2 )
  }
  cn = function(x){
    return( -(2 * A(x) - B(x) * h) )
  }
  fn = function(x){
    return( 2 * F(x) * h^2 )
  }
  
  x = seq(from = 0.1, to = 2.5, length.out = N+2)
  u = vector("numeric", N+2)
  a = sapply(x, an)
  b = sapply(x, bn)
  c = sapply(x, cn)
  f = sapply(x, fn)
  
  
  a[1] = 0
  b[1] = 1
  c[1] = 0
  f[1] = 0.8
  u[1] = 0.8
  
  a[N+2] = 0
  b[N+2] = 1
  c[N+2] = 0
  f[N+2] = 1
  u[N+2] = 1
 
  
  alpha = vector("numeric", N+2)
  beta  = vector("numeric", N+2)
  
  alpha[1] = -c[1]/b[1]
  beta[1]  = (f[1])/b[1]
  
  for(n in 2:(N + 1)){
    alpha[n]= -c[n] / (a[n] * alpha[n-1] + b[n])
    beta[n] = (f[n] - a[n] * beta[n-1]) / (a[n] * alpha[n-1] + b[n])
  }
  
  for(n in (N+1):2){
    u[n] = alpha[n] * u[n+1] + beta[n]
  }
  
  df = data.frame(x = x, y = u)
  g = ggplot(df, aes(x)) + geom_line(aes(y=y), col = "red") 
  print(g)
  print(u[N+2])
}
run(100)

вот на python:

 import math as m
import matplotlib.pyplot as plt
def run(n):
    h = 2.4 / (n+1)
    a = list(range(n)); b = list(range(n)); c = list(range(n)); f = list(range(n))
    # начальные условия
    a[0] = 0; b[0] = 1; c[0] = 0 ; f[0] = 0.8
    a[-1] = 0; b[-1] = 1; c[-1] = 0; f[-1] = 1
    alpha = list(range(n)); beta = list(range(n))
    alpha[0] = (-1)*(c[0])/(b[0])
    beta[0] = (f[0])/(b[0])
    x = list(range(n))
    x[0] = 0.1
    for i in range(1,n):
        x[i] = x[0] + i*h
    A = list(range(n)); B = list(range(n)); C = list(range(n)); F = list(range(n))
    # функци A(x), B(x), C(x), f(x).
    for i in range(n):
        A[i] = m.log(3 + m.sqrt(1+x[i]))
        B[i] = m.exp(-x[i]/2)
        C[i] = m.sqrt((3-x[i])/2)
        F[i] = 2 - m.fabs(1-x[i])
    # коэф-ты. 
    for i in range(1,n):
        a[i] = (-1) * (2 * A[i] + h * B[i])
        b[i] = (4 * A[i] + 2 * pow(h, 2) * C[i])
        c[i] = (-1) * (2 * A[i] - h * B[i])
        f[i] = 2*pow(h, 2) * f[i]
        alpha[i] = (-1)*(c[i])/(a[i]*alpha[i-1]+b[i])
        # после 5 итераций Beta[i] начинает сильно расти, а на R наооборот падает, при этом остальные 
        #                                                                          коэф-ты получаются олинаковые
        beta[i] = (f[i] - a[i]*beta[i-1])/ (a[i]*alpha[i-1]+b[i])
    # Решение u_(N+1)
    y = list(range(n))
    y[0] = 0.8
    y[-1] = 1 # кладу последний n-й элемент равный 1
    for i in reversed(range(1, n-1)):
        y[i] = alpha[i]*y[i+1]+beta[i]
    plt.plot(x, y, color='red', marker='o', linestyle='solid')
    plt.show()
run(100)

Отредактировано shevzoom (Ноя. 5, 2020 22:46:47)

py.user.next · Ноя. 6, 2020 00:04:09

Надо переводить точнее. Функция в функцию, цикл в цикл, массив в массив.

Пример

R

an = function(x) {
    return( -(2 * A(x) + B(x) * h) )
}

Python

def an(x):
    return (-(2 * A(x) + B(x) * h))

R

a = sapply(x, an)

Python

a = map(an, x)

R

for(n in (N+1):2) {
    u[n] = alpha[n] * u[n+1] + beta[n]
}

Python

for n in range(N + 1, 2 + 1):
    u[n] = alpha[n] * u[n+1] + beta[n]

И когда переводишь, смотри, что выведет на экран R в такой-то строке (вставляешь эти выводы сам в виде отладочных) и смотри что выведет Python в аналогичной строке (вставляешь эти выводы сам в виде отладочных). Когда код будет готов, отладочные выводы уберёшь.

shevzoom · Ноя. 6, 2020 18:15:27

py.user.next
Надо переводить точнее. Функция в функцию, цикл в цикл, массив в массив.ПримерR

Думаю не получится подобное в python. Я изначально делал с map() и def A(x). Будет огромное кол-во ошибок из-за этого.
Можно начать по шагам:
TypeError: ‘map’ object does not support item assignment

 import math as m
import matplotlib.pyplot as plt
def A(x):
  return( m.log(3 + m.sqrt(1 + x) ))
def B(x):
  return( m.exp(-x / 2) )
def C(x):
  return ( m.sqrt((3 - x) / 2) )
def F(x):
  return ( 2 - abs(1 - x) )
def run(n):
    h = 2 / (n+1)
    def a(x):
        return (-(2 * A(x) + B(x) * h))
    def b(x):
        return (4 * A(x) + 2 * C(x) * h ^ 2)
    def c(x):
        return (-(2 * A(x) - B(x) * h))
    def f(x):
        return (2 * F(x) * h ^ 2)
    x = list(range(n))
    y = list(range(n))
    a = map(a, x)
    b = map(a, x)
    c = map(a, x)
    f = map(a, x)
    a[0] = 0; b[0] = 1; c[0] = 0 ; f[0] = 0.4

shevzoom · Ноя. 6, 2020 18:17:42

py.user.next
Надо переводить точнее. Функция в функцию, цикл в цикл, массив в массив.ПримерR

Тоже касается и определения х, ошибка - range() не поддерживает float, можно заменить np.arrange, но тогда ломаются функции def A(x) и т.д.
Потом также циклы с np.arange() не будут считать элементы верно на каждой итерации, получаются одни и те же цифры.

py.user.next · Ноя. 6, 2020 22:40:49

shevzoom

  
TypeError: 'map' object does not support item assignment

Было так

  
a = map(a, x)

Стало так

  
a = list(map(a, x))

shevzoom
range() не поддерживает float

Делаешь функцию range_float(), которую сначала пишешь на базе своего кода, потом можешь и numpy задействовать внутри неё.

  
>>> def range_float(start, end, step, eps=0.000001):
...     x = start
...     while x <= end + eps:
...         yield x
...         x += step
... 
>>> list(range_float(2.5, 3.7, 0.2))
[2.5, 2.7, 2.9000000000000004, 3.1000000000000005, 3.3000000000000007, 3.500000000000001, 3.700000000000001]
>>>

doza_and · Ноя. 7, 2020 09:27:50

shevzoom
Будет огромное кол-во ошибок из-за этого.

Ну и что??? Если неправильно писать то в любом языке будет куча ошибок. И то что изначально написали содержит кучу ошибок.

shevzoom
Можно начать по шагам:

Это похоже правильная мысль. Но тогда ей надо следовать приводите кусок кода для одной конкретной проблемы и пишите что не так. Лично мне вообще не понятно в чем проблемы (но я не знаю R вообще).

1.

shevzoom
TypeError: ‘map’ object does not support item assignment

А зачем вам map,sapply и прочая лабуда?? Функции и так вроде векторизуются отлично в питоне.
Те проблема в том что вы используете то что вообще не нужно?

 import numpy as np
x = np.arange(0.1,10,0.5)
def A(x):
    return np.log(3+ np.sqrt(1+x))
>>> A(x)
array([1.39842273, 1.45042133, 1.4927103 , 1.52875951, 1.56039087,
       1.58869763, 1.61439474, 1.63797893, 1.6598112 , 1.68016324,
       1.69924531, 1.71722386, 1.7342332 , 1.75038346, 1.7657662 ,
       1.78045846, 1.79452574, 1.80802426, 1.82100264, 1.8335033 ])

те вместо a = sapply(x,A) просто a=A(x)

shevzoom
Потом также циклы с np.arange() не будут считать элементы верно на каждой итерации, получаются одни и те же цифры.

2.
Это вообще мне непонятно. Что вы в range задали то и будет. А как еще может быть???
Приведите пример того как вы используете range. что получается и что должно быть.

Опять. Вы используете циклы, а они тут вообще не нужны.

Питоновские циклы в numpy лучше не использовать. Их вполне заменят срезы или правильные индексы которые определены до итераций.

x = np.arange(0,6,0.5)
y = np.ones(len(x))

n = np.arange(5,2,-1)
m = np.arange(6,3,-1)
y[m] = x[n] +  y[m]
>>> y
array([1. , 1. , 1. , 1. , 2.5, 3. , 3.5, 1. , 1. , 1. , 1. , 1. ])

Отредактировано doza_and (Ноя. 7, 2020 09:38:40)

py.user.next · Ноя. 7, 2020 11:19:29

На numpy лучше не полагаться. Там может всё идти хорошо, пока не встретится какое-нибудь ограничение в numpy, которое будет невозможно обойти. Поэтому лучше написать какой-то код общего вида, который, может быть, основать на numpy потом. Но при этом, если появляется какое-то ограничение в numpy, то этот код общего вида просто отцепляется от numpy и его кусок реализуется каким-то другим образом. В итоге, для программы, использующей код в целом, будет вообще не заметно, numpy там или не numpy. Главное, что код будет выдавать всё точно на каждом шаге.

shevzoom · Ноя. 7, 2020 16:28:25

py.user.next

Действительно такая замена конструкций помогла, спасибо.
Делал отладку, писав print() к alpha, beta, y
Значения совпадали с R.
но при построение графиков на промежутке от 2 до 2.5 в R график убывает, а в python принимает значение единица, как- будто с x=2 до x=2.5 значения y(x) = 1.
Где я мог не доглядеть?

 import math as m
import matplotlib.pyplot as plt
import numpy as np
def A(x):
  return( m.log(3 + m.sqrt(1 + x) ))
def B(x):
  return( m.exp(-x / 2) )
def C(x):
  return ( m.sqrt((3 - x) / 2) )
def F(x):
  return ( 2 - abs(1 - x) )
#
# def range_float(start, end, step, eps=0.000001):
#     x = start
#     while x <= end + eps:
#         yield x
#         x += step
def run(n):
    h = 2 / (n+1)
    def an(x):
        return (-(2 * A(x) + B(x) * h))
    def bn(x):
        return (4 * A(x) + 2 * C(x) * pow(h,2))
    def cn(x):
        return (-(2 * A(x) - B(x) * h))
    def fn(x):
        return (2 * F(x) *  pow(h,2))
    x = np.arange(0.1, 2.5, h)
    a = list(map(an, x))
    b = list(map(bn, x))
    c = list(map(cn, x))
    f = list(map(fn, x))
    a[0] = 0; b[0] = 1; c[0] = 0 ; f[0] = 0.8
    a[-1] = 0; b[-1] = 1; c[-1] = 0; f[-1] = 1
    alpha = np.ones(len(x)); beta = np.ones(len(x))
    alpha[0] = -c[0]/b[0]
    beta[0] = f[0]/b[0]
    # i = 0
    for i in range(1,n):
        alpha[i] = (-1) * (c[i]) / (a[i] * alpha[i - 1] + b[i])
        beta[i] = (f[i] - a[i] * beta[i - 1]) / (a[i] * alpha[i - 1] + b[i])
    y = np.ones(len(x))
    y[0] = 0.8
    y[-1] = 1
    for i in reversed(range(1, n)):
        y[i] = alpha[i] * y[i+1] + beta[i]
        print(y[i])
    plt.plot(x, y, color='red', marker='o', linestyle='solid')
    plt.show()
run(100)

Отредактировано shevzoom (Ноя. 7, 2020 16:33:04)

py.user.next · Ноя. 7, 2020 22:49:29

shevzoom
Где я мог не доглядеть?

Например

h = 2.4 / (N + 1)

И у тебя также границы N не соответствуют в кодах на R и на Python.

Я поправлял, вывел график, как в коде на R.

  
>>> import math as m
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> 
>>> def A(x):
...   return( m.log(3 + m.sqrt(1 + x) ))
... 
>>> def B(x):
...   return( m.exp(-x / 2) )
... 
>>> def C(x):
...   return ( m.sqrt((3 - x) / 2) )
... 
>>> def F(x):
...   return ( 2 - abs(1 - x) )
... 
>>> 
>>> def range_float(start, end, step, eps=0.000001):
...     x = start
...     while x <= end + eps:
...         yield x
...         x += step
... 
>>> 
>>> def run(n):
...     h = 2.4 / (n+1)
...     def an(x):
...         return (-(2 * A(x) + B(x) * h))
...     def bn(x):
...         return (4 * A(x) + 2 * C(x) * pow(h,2))
...     def cn(x):
...         return (-(2 * A(x) - B(x) * h))
...     def fn(x):
...         return (2 * F(x) *  pow(h,2))
...     
...     #x = np.arange(0.1, 2.5, h)
...     x = list(range_float(0.1, 2.5, h))
...     print(x, len(x), n)
...     
...     a = list(map(an, x))
...     b = list(map(bn, x))
...     c = list(map(cn, x))
...     f = list(map(fn, x))
...     a[0] = 0; b[0] = 1; c[0] = 0 ; f[0] = 0.8
...     a[n+2-1] = 0; b[n+2-1] = 1; c[n+2-1] = 0; f[n+2-1] = 1
...     alpha = np.ones(len(x)); beta = np.ones(len(x))
...     alpha[0] = -c[0]/b[0]
...     beta[0] = f[0]/b[0]
...     # i = 0
...     for i in range(n + 1):
...         alpha[i] = (-1) * (c[i]) / (a[i] * alpha[i - 1] + b[i])
...         beta[i] = (f[i] - a[i] * beta[i - 1]) / (a[i] * alpha[i - 1] + b[i])
...     y = np.ones(len(x))
...     y[0] = 0.8
...     y[-1] = 1
...     for i in reversed(range(n + 1)):
...         y[i] = alpha[i] * y[i+1] + beta[i]
...         print(i, y[i])
...     plt.plot(x, y, color='red', marker='o', linestyle='solid')
...     plt.show()
... 
>>> run(100)
[0.1, 0.12376237623762376, 0.14752475247524752, 0.17128712871287127, 0.19504950495049503, 0.21881188118811878, 0.24257425742574254, 0.2663366336633663, 0.2900990099009901, 0.31386138613861386, 0.33762376237623765, 0.36138613861386143, 0.3851485148514852, 0.408910891089109, 0.4326732673267328, 0.45643564356435656, 0.48019801980198035, 0.5039603960396041, 0.5277227722772279, 0.5514851485148516, 0.5752475247524754, 0.5990099009900992, 0.622772277227723, 0.6465346534653468, 0.6702970297029706, 0.6940594059405943, 0.7178217821782181, 0.7415841584158419, 0.7653465346534657, 0.7891089108910895, 0.8128712871287133, 0.836633663366337, 0.8603960396039608, 0.8841584158415846, 0.9079207920792084, 0.9316831683168322, 0.955445544554456, 0.9792079207920797, 1.0029702970297034, 1.0267326732673272, 1.050495049504951, 1.0742574257425748, 1.0980198019801986, 1.1217821782178223, 1.1455445544554461, 1.16930693069307, 1.1930693069306937, 1.2168316831683175, 1.2405940594059413, 1.264356435643565, 1.2881188118811888, 1.3118811881188126, 1.3356435643564364, 1.3594059405940602, 1.383168316831684, 1.4069306930693077, 1.4306930693069315, 1.4544554455445553, 1.478217821782179, 1.5019801980198029, 1.5257425742574267, 1.5495049504950504, 1.5732673267326742, 1.597029702970298, 1.6207920792079218, 1.6445544554455456, 1.6683168316831694, 1.6920792079207931, 1.715841584158417, 1.7396039603960407, 1.7633663366336645, 1.7871287128712883, 1.810891089108912, 1.8346534653465358, 1.8584158415841596, 1.8821782178217834, 1.9059405940594072, 1.929702970297031, 1.9534653465346548, 1.9772277227722785, 2.0009900990099023, 2.024752475247526, 2.04851485148515, 2.0722772277227737, 2.0960396039603975, 2.1198019801980212, 2.143564356435645, 2.167326732673269, 2.1910891089108926, 2.2148514851485164, 2.23861386138614, 2.262376237623764, 2.2861386138613877, 2.3099009900990115, 2.3336633663366353, 2.357425742574259, 2.381188118811883, 2.4049504950495066, 2.4287128712871304, 2.452475247524754, 2.476237623762378, 2.5000000000000018] 102 100
100 1.00478041936
99 1.00953665589
98 1.01426507668
97 1.01896200308
96 1.02362371461
95 1.02824645232
94 1.03282642181
93 1.03735979586
92 1.04184271668
91 1.04627129793
90 1.05064162649
89 1.05494976396
88 1.05919174799
87 1.06336359342
86 1.06746129324
85 1.0714808194
84 1.07541812352
83 1.07926913745
82 1.08302977369
81 1.08669592576
80 1.09026346844
79 1.09372825795
78 1.09708613203
77 1.10033290992
76 1.10346439238
75 1.10647636149
74 1.10936458054
73 1.11212479381
72 1.11475272624
71 1.11724408318
70 1.11959454999
69 1.12179979169
68 1.1238554525
67 1.12575715539
66 1.12750050162
65 1.1290810702
64 1.13049441738
63 1.13173607609
62 1.13280155536
61 1.13368633975
60 1.13438588873
59 1.13489563607
58 1.13521098923
57 1.13532732873
56 1.13524000747
55 1.13494435014
54 1.13443565254
53 1.13370918094
52 1.13276017141
51 1.13158382921
50 1.13017532814
49 1.12852980986
48 1.12664238329
47 1.12450812398
46 1.12212207349
45 1.11947923878
44 1.11657459161
43 1.11340306798
42 1.10995956753
41 1.10623895303
40 1.1022360498
39 1.09794564524
38 1.09336248831
37 1.08848128905
36 1.08331247225
35 1.07786856524
34 1.07216192969
33 1.06620476327
32 1.06000910119
31 1.05358681794
30 1.04694962893
29 1.04010909222
28 1.03307661028
27 1.02586343174
26 1.01848065318
25 1.01093922093
24 1.00324993288
23 0.995423440326
22 0.987470249781
21 0.979400724824
20 0.971225087952
19 0.962953422419
18 0.954595674092
17 0.946161653285
16 0.937661036609
15 0.929103368804
14 0.92049806456
13 0.911854410343
12 0.90318156619
11 0.89448856751
10 0.885784326854
9 0.87707763568
8 0.868377166095
7 0.85969147258
6 0.851028993693
5 0.842398053754
4 0.833806864505
3 0.825263526745
2 0.816776031948
1 0.808352263847
0 0.8
>>>

shevzoom · Ноя. 8, 2020 13:20:27

py.user.next

Видимо, не заметил на стер “.4” , моя невнимательность, спасибо за помощь.

Python-сообщество

Уведомления

#1 Ноя. 5, 2020 22:44:35

Из-за чего происходит разница в графиках по оси ординат в R и Python?

#2 Ноя. 6, 2020 00:04:09

Из-за чего происходит разница в графиках по оси ординат в R и Python?

#3 Ноя. 6, 2020 18:15:27

Из-за чего происходит разница в графиках по оси ординат в R и Python?

#4 Ноя. 6, 2020 18:17:42

Из-за чего происходит разница в графиках по оси ординат в R и Python?

#5 Ноя. 6, 2020 22:40:49

Из-за чего происходит разница в графиках по оси ординат в R и Python?

#6 Ноя. 7, 2020 09:27:50

Из-за чего происходит разница в графиках по оси ординат в R и Python?

#7 Ноя. 7, 2020 11:19:29

Из-за чего происходит разница в графиках по оси ординат в R и Python?

#8 Ноя. 7, 2020 16:28:25

Из-за чего происходит разница в графиках по оси ординат в R и Python?

#9 Ноя. 7, 2020 22:49:29

Из-за чего происходит разница в графиках по оси ординат в R и Python?

#10 Ноя. 8, 2020 13:20:27

Из-за чего происходит разница в графиках по оси ординат в R и Python?

Board footer