Почему Matlab и Python хранят значения float64 по-разному? - PullRequest
Новая
       16

Почему Matlab и Python хранят значения float64 по-разному?

0 голосов
/ 01 ноября 2019

В настоящее время я решаю двумерный ODE с использованием Python. В прошлом году я написал тот же алгоритм в Matlab для своего проекта, и он дал мне желаемые результаты. Однако, хотя алгоритм тот же, Python, похоже, хранит значения float64 немного по-разному (т. Е. Разные десятичные значения, например, третий элемент массива N1 в Matlab - 3009.570691672611, а в Python - 3009.5706916726026), и после всех итераций он идетчерез это в конечном итоге дает совсем другие результаты, чем код Matlab. Почему это, и как мне это исправить?

Вот код Matlab: 

function y = treatiseOCP
format long
test = -1;
T = 5;
delta = 0.001;
M = 999;
t = linspace(0,T,M+1);
h = T/M;
h2 = h/2;
C = 100;
uMaximum = 0.8;
%Parameter values

r=0.84;
K=5000;
g=0.2833;
c=0.6;
a=2.5;
b=0.8;
d=0.7;
m=0.0625;

N1=zeros(1,M+1);
N2=zeros(1,M+1);

u = zeros(1,M+1);

lambda1 = zeros(1,M+1);
lambda1(M+1) = -1;

lambda2 = zeros(1,M+1);
lambda2(M+1) = 1;

N1(1) = 3000;
N2(1) = 300;

while(test < 0)

    oldu = u;
    oldN1 = N1;
    oldN2 = N2;
    oldlambda1 = lambda1;
    oldlambda2 = lambda2;

    for i = 1:M
    %k1,k2,k3,k4 values for N1 and N2 respectively
        k11 = r*N1(i)*(1-N1(i)/K) - (g*c*N1(i)*N2(i))/(a+N1(i));
        k12 = (b*d*N1(i)*N2(i))/(a+N1(i)) - m*N2(i) - u(i)*N2(i); 

        k21 = r*(N1(i)+h2*k11)*(1-(N1(i)+h2*k11)/K) - (g*c*(N1(i)+h2*k11)*(N2(i)+h2*k12))/...
             (a+(N1(i)+h2*k11));
        k22 = (b*d*(N1(i)+h2*k11)*(N2(i)+h2*k12))/(a+(N1(i)+h2*k11))...
             - m*(N2(i)+h2*k12) - 0.5*(u(i)+u(i+1))*(N2(i)+h2*k12);


        k31 = r*(N1(i)+h2*k21)*(1-(N1(i)+h2*k21)/K) - (g*c*(N1(i)+h2*k21)*(N2(i)+h2*k22))/...
             (a+(N1(i)+h2*k21));
        k32 = (b*d*(N1(i)+h2*k21)*(N2(i)+h2*k22))/(a+(N1(i)+h2*k21))...
             - m*(N2(i)+h2*k22) - 0.5*(u(i)+u(i+1))*(N2(i)+h2*k22);

        k41 = r*(N1(i)+h*k31)*(1-(N1(i)+h*k31)/K) - (g*c*(N1(i)+h*k31)*(N2(i)+h*k32))/...
             (a+(N1(i)+h*k31));
        k42 = (b*d*(N1(i)+h*k31)*(N2(i)+h*k32))/(a+(N1(i)+h*k31))...
             - m*(N2(i)+h*k32) - u(i+1)*(N2(i)+h*k32);

        N1(i+1) = N1(i) + (h/6)*(k11 + 2*k21 + 2*k31 + k41);
        N2(i+1) = N2(i) + (h/6)*(k12 + 2*k22 + 2*k32 + k42);

    end

    for i = 1:M
    %n1,n2,n3,n4 values for lambda1 and lambda2 respectively   
        j = M + 2 - i;

        n11 = -lambda1(j)*(r-(2*r/K)*N1(j) - g*c*(((a+N1(j))*N2(j) - N1(j)*N2(j))/(a+N1(j))^2))...
              -lambda2(j)*b*d*(((a+N1(j))*N2(j) - N1(j)*N2(j))/(a+N1(j))^2);
        n12 = -lambda1(j)*(-g*c*N1(j)/(a+N1(j))) - lambda2(j)*(b*d*N1(j)/(a+N1(j)) - m - u(j));

        n21 = -(lambda1(j)-h2*n11)*(r-(2*r/K)*0.5*(N1(j)+N1(j-1)) - g*c*(((a+0.5*(N1(j)+N1(j-1)))...
              *0.5*(N2(j)+N2(j-1)) - 0.5*(N1(j)+N1(j-1))*0.5*(N2(j)+N2(j-1))))/...
              (a+0.5*(N1(j)+N1(j-1)))^2)...
              -(lambda2(j)-h2*n12)*b*d*(((a+0.5*(N1(j)+N1(j-1)))*0.5*(N2(j)+N2(j-1))...
              -0.5*(N1(j)+N1(j-1))*0.5*(N2(j)+N2(j-1)))/(a+0.5*(N1(j)+N1(j-1)))^2);
        n22 =  -(lambda1(j)+h2*n11)*(-g*c*0.5*(N1(j)+N1(j-1))/(a+0.5*(N1(j)+N1(j-1))))...
               - (lambda2(j)+h2*n12)*(b*d*0.5*(N1(j)+N1(j-1))/(a+0.5*(N1(j)+N1(j-1))) - m - u(j));

        n31 = -(lambda1(j)-h2*n21)*(r-(2*r/K)*0.5*(N1(j)+N1(j-1)) - g*c*(((a+0.5*(N1(j)+N1(j-1)))...
              *0.5*(N2(j)+N2(j-1)) - 0.5*(N1(j)+N1(j-1))*0.5*(N2(j)+N2(j-1))))/...
              (a+0.5*(N1(j)+N1(j-1)))^2)...
              -(lambda2(j)-h2*n22)*b*d*(((a+0.5*(N1(j)+N1(j-1)))*0.5*(N2(j)+N2(j-1))...
              -0.5*(N1(j)+N1(j-1))*0.5*(N2(j)+N2(j-1)))/(a+0.5*(N1(j)+N1(j-1)))^2);
        n32 =  -(lambda1(j)+h2*n21)*(-g*c*0.5*(N1(j)+N1(j-1))/(a+0.5*(N1(j)+N1(j-1))))...
               - (lambda2(j)+h2*n22)*(b*d*0.5*(N1(j)+N1(j-1))/(a+0.5*(N1(j)+N1(j-1))) - m - u(j));

        n41 = -(lambda1(j)-h*n31)*(r-(2*r/K)*N1(j-1) - g*c*(((a+N1(j-1))...
              *N2(j-1) - N1(j-1)*N2(j-1)))/...
              (a+N1(j-1))^2)...
              -(lambda2(j)-h*n32)*b*d*(((a+N1(j-1))*N2(j-1)...
              -N1(j-1)*N2(j-1))/(a+N1(j-1))^2);
        n42 = -(lambda1(j)+h*n31)*(-g*c*N1(j-1)/(a+N1(j-1)))...
               - (lambda2(j)+h*n32)*(b*d*N1(j-1)/(a+N1(j-1)) - m - u(j));


        lambda1(j-1) = lambda1(j) - h/6*(n11 + 2*n21 + 2*n31 + n41);
        lambda2(j-1) = lambda2(j) - h/6*(n12 + 2*n22 + 2*n32 + n42);

    end

    u1 = min(uMaximum,max(0,lambda2.*N2/C));
    u = 0.5*(u1 + oldu);


    temp1 = delta*sum(abs(u)) - sum(abs(oldu - u));
    temp2 = delta*sum(abs(N1)) - sum(abs(oldN1 - N1));
    temp3 = delta*sum(abs(N2)) - sum(abs(oldN2 - N2));
    temp4 = delta*sum(abs(lambda1)) - sum(abs(oldlambda1 - lambda1));
    temp5 = delta*sum(abs(lambda2)) - sum(abs(oldlambda2 - lambda2));

    test = min(temp1, min(temp2, min(temp3, min(temp4, min(temp5, min(temp4, temp5))))));
end

y(1,:) = t;
y(2,:) = N1;
y(3,:) = N2;
y(4,:) = lambda1;
y(5,:) = lambda2;
y(6,:) = u;

Z=[y(1,:); y(2, :); y(3, :); y(4, :); y(5, :); y(6, :)];
fid = fopen('myresults_TreatiseOCP.txt','w');
fprintf(fid,'%f %f %f %f %f %f\r\n',Z);
fclose(fid);

figure(1)
            subplot(4,1,1);plot(y(1,:),y(2,:),'LineWidth',2)
            subplot(4,1,1);xlabel('Time (years)')
            subplot(4,1,1);ylabel('Prey')
            subplot(4,1,2);plot(y(1,:),y(3,:),'LineWidth',2)
            subplot(4,1,2);xlabel('Time (years)')
            subplot(4,1,2);ylabel('Predator')
            subplot(4,1,3);plot(y(1,:),y(6,:),'LineWidth',2)
            subplot(4,1,3);xlabel('Time (years)')
            subplot(4,1,3);ylabel('Harvesting rate')

А вот код Python: 


import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

test = -1
T = 5 #years
delta = 0.001
M = 999
t = np.linspace(0,T,num=M+1).reshape((1,M+1))
h = T/M #step-size
h2 = h/2
C = 100
uMaximum = 0.8

#Parameter values
r=0.84
K=5000
g=0.2833
c=0.6
a=2.5
b=0.8
d=0.7
m=0.0625

N1 = np.zeros(M+1).reshape((1,M+1))
N2 = np.zeros(M+1).reshape((1,M+1))
u = np.zeros(M+1).reshape((1,M+1))

lambda1 = np.zeros(M+1).reshape((1,M+1))
lambda1[:,-1] = -1
lambda2 = np.zeros(M+1).reshape((1,M+1))
lambda2[:,-1] = 1

N1[:,0] = 3000
N2[:,0] = 300
#all above gives same as Matlab code

while test < 0:
    oldu = u.copy()
    oldN1 = N1.copy()
    oldN2 = N2.copy()
    oldlambda1 = lambda1.copy()
    oldlambda2 = lambda2.copy()

    for i in range(M):
        # k1,k2,k3,k4 values for N1 and N2 respectively
        k11 = r * N1[0,i] * (1 - N1[0,i] / K) - (g * c * N1[0,i] * N2[0,i]) / (a + N1[0,i])
        k12 = (b * d * N1[0,i] * N2[0,i]) / (a + N1[0,i]) - m * N2[0,i] - u[0,i] * N2[0,i]

        k21 = r * (N1[0,i] + h2 * k11) * (1 - (N1[0,i] + h2 * k11) / K) - (g * c * (N1[0,i] + h2 * k11) * (N2[0,i] + h2 * k12))\
              /(a + (N1[0,i] + h2 * k11))
        k22 = (b * d * (N1[0,i] + h2 * k11) * (N2[0,i] + h2 * k12)) / (a + (N1[0,i] + h2 * k11)) \
              - m * (N2[0,i] + h2 * k12) - 0.5 * (u[0,i] + u[0,i+1]) * (N2[0,i] + h2 * k12)

        k31 = r * (N1[0,i] + h2 * k21) * (1 - (N1[0,i] + h2 * k21) / K) - (g * c * (N1[0,i] + h2 * k21) * (N2[0,i] + h2 * k22))\
              /(a + (N1[0,i] + h2 * k21))
        k32 = (b * d * (N1[0,i] + h2 * k21) * (N2[0,i] + h2 * k22)) / (a + (N1[0,i] + h2 * k21)) \
              - m * (N2[0,i] + h2 * k22) - 0.5 * (u[0,i] + u[0,i+1]) * (N2[0,i] + h2 * k22)

        k41 = r * (N1[0,i] + h * k31) * (1 - (N1[0,i] + h * k31) / K) - (g * c * (N1[0,i] + h * k31) * (N2[0,i] + h * k32)) / \
              (a + (N1[0,i] + h * k31))
        k42 = (b * d * (N1[0,i] + h * k31) * (N2[0,i] + h * k32)) / (a + (N1[0,i] + h * k31)) \
              - m * (N2[0,i] + h * k32) - u[0,i+1] * (N2[0,i] + h * k32)

        N1[0,i + 1] = N1[0,i] + (h / 6) * (k11 + 2 * k21 + 2 * k31 + k41)
        N2[0,i + 1] = N2[0,i] + (h / 6) * (k12 + 2 * k22 + 2 * k32 + k42)


    for i in range(M):
        # n1,n2,n3,n4 values for lambda1 and lambda2 respectively
        j = M-i
        n11 = -lambda1[0,j] * (r - (2 * r / K) * N1[0,j] - g * c * ((a + N1[0,j]) * N1[0,j] - N1[0,j] * N1[0,j])\
                /(a + np.power(N1[0,j], 2))) - lambda2[0,j] * b * d * ((a + N1[0,j]) * N1[0,j] - N1[0,j] * N1[0,j])\
              / np.power(a + N1[0,j], 2)
        n12 = -lambda1[0,j] * (-g * c * N1[0,j] / (a + N1[0,j]) - lambda2[0,j] * (b * d * N1[0,j] / (a + N1[0,j]) - m - u[0,j]))

        n21 = -(lambda1[0,j] - h2 * n11) * (r - (2 * r / K) * 0.5 * (N1[0,j] + N1[0,j-1]) - g * c *\
                (((a + 0.5 * (N1[0,j] + N1[0,j-1])) * 0.5 * (N1[0,j] + N2[0,j-1]) - 0.5 * (N1[0,j] + N1[0,j-1])\
                  * 0.5 * (N1[0,j] + N2[0,j-1]))) / np.power(a + 0.5 * (N1[0,j] + N1[0,j-1]), 2)) \
             - (lambda2[0,j] - h2 * n12) * b * d * (((a + 0.5 * (N1[0,j] + N1[0,j-1])) * 0.5 * (N1[0,j] + N2[0,j-1]) \
               - 0.5 * (N1[0,j] + N1[0,j-1]) * 0.5 * (N1[0,j] + N2[0,j-1])) / np.power(a + 0.5 * (N1[0,j] + N1[0,j-1]), 2))
        n22 = -(lambda1[0,j]+h2 * n11) * (-g * c * 0.5 * (N1[0,j] + N1[0,j-1]) / (a + 0.5 * (N1[0,j] + N1[0,j-1])))\
                - (lambda2[0,j] + h2 * n12) * (b * d * 0.5 * (N1[0,j] + N1[0,j-1]) / (a + 0.5 * (N1[0,j] + N1[0,j-1])) - m - u[0,j])

        n31 = -(lambda1[0,j]-h2*n21)*(r-(2*r/K)*0.5*(N1[0,j]+N1[0,j-1]) - g*c*(((a+0.5*(N1[0,j]+N1[0,j-1]))*0.5*(N1[0,j]+N2[0,j-1])\
                - 0.5*(N1[0,j]+N1[0,j-1])*0.5*(N1[0,j]+N2[0,j-1])))/np.power(a+0.5*(N1[0,j]+N1[0,j-1]),2))\
                -(lambda2[0,j]-h2*n22)*b*d*(((a+0.5*(N1[0,j]+N1[0,j-1]))*0.5*(N1[0,j]+N2[0,j-1])\
                -0.5*(N1[0,j]+N1[0,j-1])*0.5*(N1[0,j]+N2[0,j-1]))/np.power(a+0.5*(N1[0,j]+N1[0,j-1]),2))
        n32 =  -(lambda1[0,j]+h2 * n21) * (-g * c * 0.5 * (N1[0,j] + N1[0,j-1]) / (a + 0.5 * (N1[0,j] + N1[0,j-1]))) \
                - (lambda2[0,j] + h2 * n22) * (b * d * 0.5 * (N1[0,j] + N1[0,j-1]) / (a + 0.5 * (N1[0,j] + N1[0,j-1])) - m - u[0,j])

        n41 = -(lambda1[0,j] - h * n31) * (r - (2 * r / K) * N1[0,j-1] - g * c * (((a + N1[0,j-1]) \
                * N2[0,j-1] - N1[0,j-1] * N2[0,j-1])) / np.power(a + N1[0,j-1], 2)) \
              - (lambda2[0,j] - h * n32) * b * d * (((a + N1[0,j-1]) * N2[0,j-1] \
                - N1[0,j-1] * N2[0,j-1]) / np.power(a + N1[0,j-1], 2))
        n42 = -(lambda1[0,j] + h * n31) * (-g * c * N1[0,j-1] / (a + N1[0,j-1])) \
              - (lambda2[0,j] + h * n32) * (b * d * N1[0,j-1] / (a + N1[0,j-1]) - m - u[0,j])

        lambda1[0,j - 1] = lambda1[0,j] - h / 6 * (n11 + 2 * n21 + 2 * n31 + n41)
        lambda2[0,j - 1] = lambda2[0,j] - h / 6 * (n12 + 2 * n22 + 2 * n32 + n42)
        u1 = min(uMaximum, max(0, np.max(lambda2 * N2 / C)))

        u = 0.5 * (u1 + oldu)

        temp1 = delta * np.sum(abs(u)) - np.sum(abs(oldu - u))
        temp2 = delta * np.sum(abs(N1)) - np.sum(abs(oldN1 - N1))
        temp3 = delta * np.sum(abs(N2)) - np.sum(abs(oldN2 - N2))
        temp4 = delta * np.sum(abs(lambda1)) - np.sum(abs(oldlambda1 - lambda1))
        temp5 = delta * np.sum(abs(lambda2)) - np.sum(abs(oldlambda2 - lambda2))

        test = min(temp1, min(temp2, min(temp3, min(temp4, min(temp5, min(temp4, temp5))))))
y = np.zeros((6,M+1))
y[0,:] = t
y[1,:] = N1
y[2,:] = N2
y[3,:] = lambda1
y[4,:] = lambda2
y[5,:] = u

Z = np.matrix(y)
with open('myresults_Python.txt','wb') as f:
    for line in Z:
        np.savetxt(f, line, fmt='%.6f')

y = pd.DataFrame({'t': y[0,:], 'N1': y[1,:], 'N2' : y[2,:], 'lamda1' : y[3,:],\
                  'lambda2' : y[4,:], 'u' : y[5,:]})


sns.lineplot(x='t',y='N1',data=y)
plt.xlabel('Time (years)')
plt.ylabel('Prey')
plt.xlim((0,5))
plt.ylim((3000,5000))

sns.lineplot(x='t',y='N2',data=y)
plt.xlabel('Time (years)')
plt.ylabel('Predator')
plt.xlim((0,5))
plt.ylim((50,300))

sns.lineplot(x='t',y='u',data=y)
plt.xlabel('Time (years)')
plt.ylabel('Harvesting rate')
plt.xlim((0,5))
plt.ylim((0.68,0.8))

Вы можетеПосмотрите на визуализации (графики), что они немного отличаются, особенно термин «Уборка урожая». Я ожидаю, что они будут одинаковыми.

Важно отметить, что некоторые значения хранятся с одинаковой точностью на обоих языках, но некоторые нет, и я хотел бы знать, почему?

...