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

shevzoom · Дек. 14, 2020 20:36:37

Добрый день, программирую уравнение переноса, и испытываю трудности с заполнением матрицы.
Где-то видимо я не доглядел в индексах, не понимаю где.
Получаю: “ index 380 is out of bounds for axis 0 with size 380 ”
это в моменте цикл в цикле U.

     U = np.zeros((n, m))
    for i in range(0, n):
        U[1, ] = np.exp(-x[i] / 2)
    for j in range(0, n-1):
        for i in range(n, m):
            U[i+1, j] = (1+c*r)*U[i,j] - c*r*U[i,j-1] + tau*(np.log(1+t[i]) - 2*x[i])

вот полный код:

 import numpy as np
import matplotlib.pyplot as plt
import math as ma
def solve():
    h = 0.001
    c = 0.4
    T = 1.9
 
    xmin = 0.2
    xmax = 2.2
    
    h = 0.01
    r = 0.5
    tau = r*h
    n = ma.floor(T/tau)
    #right bound
    m = ma.floor((xmax - xmin)/h) + n
    t = []
    
    for i in range(n): 
      t.append(i*tau)
    x = np.linspace(xmin, xmax, m)
    U = np.zeros((n, m))
    
    for i in range(0, n):
        U[1, ] = np.exp(-x[i] / 2)
        
    
    for j in range(0, n-1):
        for i in range(n, m):
            U[i+1, j] = (1+c*r)*U[i,j] - c*r*U[i,j-1] + tau*(np.log(1+t[i]) - 2*x[i])
    y = np.ones(len(x))
    for i in range(0,n):
        y[i] = pow((x[i]-2.5*t[i]),2) - t[i] + (1+t[i])*np.log(t[i]+1) + np.exp(-(x[i]-2.5*t[i])/2) - pow(x[i],2)
    plt.plot(x, U[:,i])
    plt.plot(x, y)
    #plt.xlim(xmin, xmax)
    #plt.ylim(0, T)
    plt.show()
    
solve()

вот, например, та же самая задача на R, работает. Я хочу упростить немного, переписав на python, поэтому код не 1 в 1 копировал, но по идее суть остается одна.

 require(ggplot2)
Solve = function(h, r){
    if(r > 1)
      return("��� r > 1 ����� ����������")
    if(r <= 0)
      return("r ������ ���� ������ 0")
    
    phiFunc = function(x, t){
  #return(x^2/2 + 2*t)
  return(x^2/2 + 2*t)
  exp(-x/2)
}
    psi = function(x){
  #return(1 + cos(x))
  return(exp(-x/2))
}
    u_an = function(x, t){
      #return((x^3)/6 - ((x-t)^3)/6 + t^2 + cos(x-t) + 1)
      return((x-2.5*t)^2 - t + (1+t)*log(1+t) + exp(-(x-2.5*t)/2) - (x^2))
    }
    
    # a = -1
    # b =  1
    # c = -1
    # T =  2
    a = 0.2
    b =  2.2
    c = -0.4
    T =  1.9
    
    #r = 1
    #h = 0.01
    tau = r * h
    
    #up bound
    N = floor(T/tau)
    #right bound
    M = floor((b - a)/h) + N
    
    x = seq(from = b - M*h, length.out = M, by = h)
    u = matrix(0, N, M)
    
    u[1, ] = sapply(x, psi)
    
    for(n in 1:(N - 1)){
        phi = sapply(x, function(x){ return(phiFunc(x, tau * (n-1))) })
        for(m in (n + 1):M){
            u[n + 1, m] = tau * phi[m] + u[n, m] + c*r*(u[n, m] - u[n, m-1])
        }
        if(n%%100 == 0) {
          print(sprintf("%.2f%%", 100*n/N)) 
        }
    }
    
    
    df = data.frame(x  = x[(N+1):M],
                    u1 = u[    1, (N+1):M], 
                    u2 = u[  N/3, (N+1):M],
                    u3 = u[2*N/3, (N+1):M],
                    u4 = u[  N-1, (N+1):M],
                    u_an1 = sapply(x, function(x) u_an(x,   0)  ) [(N+1):M],
                    u_an2 = sapply(x, function(x) u_an(x, T/3)  ) [(N+1):M],
                    u_an3 = sapply(x, function(x) u_an(x, 2*T/3)) [(N+1):M],
                    u_an4 = sapply(x, function(x) u_an(x,   T)  ) [(N+1):M]
    )
    
    g = ggplot(df, aes(x)) + geom_line(aes(y = u1)) + geom_line(aes(y = u_an1), col = 2) 
    print(g)
    g = ggplot(df, aes(x)) + geom_line(aes(y = u2)) + geom_line(aes(y = u_an2), col = 2) 
    print(g)
    g = ggplot(df, aes(x)) + geom_line(aes(y = u3)) + geom_line(aes(y = u_an3), col = 2) 
    print(g)
    g = ggplot(df, aes(x)) + geom_line(aes(y = u4)) + geom_line(aes(y = u_an4), col = 2) 
    print(g)
    print("100%")
}  
Solve(0.01, 0.5)

UPD: Кажется, ошибка где-то еще, моя разностная схема, не сходится к точному решению на всех временных слоях кроме нулевого, на нулевом сходится. Вроде бы пси пихаем и нормально, а потом куда-то расходится. На этом форуме, есть раздел, где можно задать вопросы по математике?

Отредактировано shevzoom (Дек. 14, 2020 20:51:07)

py.user.next · Дек. 14, 2020 23:00:15

shevzoom
Где-то видимо я не доглядел в индексах, не понимаю где.

Ну, это было и в прошлый раз, когда ты циклы неправильно переводил. Я показывал, как надо переводить for на языке R в for на языке Python. Это разные конструкции, хоть у них и одинаковое слово for. Диапазоны в языке R не равны диапазонам в языке Python. Поэтому всё это записывается по-разному.

shevzoom
На этом форуме, есть раздел, где можно задать вопросы по математике?

Есть форумы по математике. Крупный форум https://dxdy.ru/ .

shevzoom · Дек. 15, 2020 09:12:42

Да, я знаю, что в R цикл идет по индексам элементов, а в python по элементам.
хорошо переделаю как в прошлый раз, а как перевести это, не подскажете?

   for(n in 1:(N - 1)){
        phi = sapply(x, function(x){ return(phiFunc(x, tau * (n-1))) })

py.user.next

Отредактировано shevzoom (Дек. 15, 2020 09:18:28)

py.user.next · Дек. 15, 2020 09:19:10

shevzoom
Да, я знаю, что в R цикл идет по индексам элементов, а в python по элементам.

Речь вообще не об этом идёт. У диапазонов в R и диапазонов в Python разная семантика. В R диапазон 1:5 даёт последовательность 1 2 3 4 5, а в Python диапазон 1:5 даёт последовательность 1 2 3 4, потому что правая граница - это граница остановки. Так же и срезы в Python устроены.

В языке Python

  
>>> list(range(1, 5))
[1, 2, 3, 4]
>>> 
>>> list(range(5, 1))
[]
>>> 
>>> list(range(5, 1, -1))
[5, 4, 3, 2]
>>>
>>> list(range(5))
[0, 1, 2, 3, 4]
>>>
>>> [1, 2, 3, 4, 5, 6, 7, 8, 9, 10][1:5]  # правая граница останавливает выборку на элементе 6, который имеет индекс 5 в списке
[2, 3, 4, 5]
>>>

В языке R

> 1:5
[1] 1 2 3 4 5
> 
> 10:20
 [1] 10 11 12 13 14 15 16 17 18 19 20
> 
> 5:1
[1] 5 4 3 2 1
>

Как видишь, вообще по-разному идёт генерация последовательности. Поэтому и нельзя напрямую переводить буква в букву, как ты это делаешь. Надо понимать, что там генерируется, а для этого надо везде вставлять print'ы, чтобы на экране писалось содержимое каждой последовательности и каждой переменной и таким образом видеть, правильно ты перевёл или вообще не то выдаётся. Не надо думать; надо смотреть конкретно; а для этого есть оператор вывода содержимого переменной или последовательности на экран.

shevzoom
а как перевести это

for(n in 1:(N - 1)){
        phi = sapply(x, function(x){ return(phiFunc(x, tau * (n-1))) })

  
for n in range(1, (N - 1) + 1):
    phi = list(map(lambda v: phiFunc(v, tau * (n - 1)), x))

Отредактировано py.user.next (Дек. 15, 2020 09:38:31)

shevzoom · Дек. 15, 2020 21:33:00

py.user.next

Я переписывал, как мы договорились в прошлый раз, но уже получил
ValueError: cannot copy sequence with size 58361 to array axis with dimension 580.
Аналог matrix(0, N, M) ведь np.zeros((N, M)) ведь?

 import numpy as np
import matplotlib.pyplot as plt
import math as ma
def range_float(start, end, step, eps=0.000001):
    x = start
    while x <= end + eps:
        yield x
        x += step
def phiF(x,t):
    return(log(1+t) - 2*x)
def psi(x):
    return(ma.exp(-x/2))
def u_n(x,t):
    return((x-2.5*t)**2 - t + (1+t)*log(1+t) + exp(-(x-2.5*t)/2) - (x**2))
    
def solve(h, r):
    a = 0.2
    b =  2.2
    c = -0.4
    T =  1.9
 
    tau = r*h
    N = ma.floor(T/tau)
    M = ma.floor((b - a)/h) + N
    x = list(range_float(b - M*h, M, h))
    u = np.zeros((N, M))
    
    u[0, ] = list(map(psi, x))
#     u[1, ] = list(map(psi, x))
#     u[1, :] = list(map(psi, x))
   
        
    for n in range(1, (N - 1) + 1):
        phi = list(map(lambda v: phiFunc(v, tau * (n - 1)), x))
        for m in range((n+1)+1, m+1):
            u[n + 1, m] = tau * phi[m] + u[n, m] + c*r*(u[n, m] - u[n, m-1])
    plt.plot(x, U[:,i])
    plt.show()
    
solve(0.01, 0.5)

Отредактировано shevzoom (Дек. 15, 2020 21:35:53)

py.user.next · Дек. 15, 2020 22:06:30

Вот ты пишешь

shevzoom

  
for m in range((n+1)+1, m+1):

А там написано

shevzoom

for(m in (n + 1):M){

На Python оно будет выглядеть вот так

  
for m in range(n + 1, M + 1):

Вот теперь со своим сравни и подумай, в состоянии ли ты вообще перевести хоть что-нибудь.

shevzoom
Я переписывал, как мы договорились в прошлый раз, но уже получил
ValueError: cannot copy sequence with size 58361 to array axis with dimension 580.
Аналог matrix(0, N, M) ведь np.zeros((N, M)) ведь?

Это тут вообще ни при чём. Ты не в состоянии просто буквы переписать. Ты даже не понимаешь, что большие буквы и маленькие буквы - это разные буквы. Переменная m и переменная M - это разные переменные.

Так что, с чего ты взял, что ты вообще можешь этот код переписать с R на какой-то другой язык, если ты не знаешь ни R, ни какого-то другого языка?

Отредактировано py.user.next (Дек. 15, 2020 22:08:29)

shevzoom · Дек. 16, 2020 09:08:25

py.user.next
u = list(map(psi, x))

Причем тут вообще этот цикл? Он для заполнения всех строк матрицы, кроме 1-й, я даже его не смотрел еще.
Ошибка происходит раньше , а конкретнее:

     u[0, ] = list(map(psi, x))

Отредактировано shevzoom (Дек. 16, 2020 18:04:18)

py.user.next · Дек. 16, 2020 23:44:39

shevzoom
Причем тут вообще этот цикл? Он для заполнения всех строк матрицы, кроме 1-й, я даже его не смотрел еще.

Цикл неправильно переведён с R на Python. Его надо не смотреть, его надо сначала правильно переводить, как и всё остальное, а уже потом только смотреть, где остались недочёты.

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

Уведомления

#1 Дек. 14, 2020 20:36:37

Почему возникает ошибка циклов?

#2 Дек. 14, 2020 23:00:15

Почему возникает ошибка циклов?

#3 Дек. 15, 2020 09:12:42

Почему возникает ошибка циклов?

#4 Дек. 15, 2020 09:19:10

Почему возникает ошибка циклов?

#5 Дек. 15, 2020 21:33:00

Почему возникает ошибка циклов?

#6 Дек. 15, 2020 22:06:30

Почему возникает ошибка циклов?

#7 Дек. 16, 2020 09:08:25

Почему возникает ошибка циклов?

#8 Дек. 16, 2020 23:44:39

Почему возникает ошибка циклов?

Board footer