Одно и то же уравнение, разные ответы из Pylab и Octave

#python #matlab #octave #equation #matplotlib

#python #matlab #octave #уравнение #matplotlib

Вопрос:

Я переношу код, созданный в octave, в pylab. Одно из перенесенных уравнений дает совершенно иные результаты в python, чем в octave.

Лучший способ объяснить — показать графики, сгенерированные octave и pylab из одного и того же уравнения.

Вот упрощенный фрагмент исходного уравнения в octave. В этом небольшом тестовом скрипте результат функции с нулевым значением phi выводится из ~ (-pi, pi):

 clear
clc
close all

L1 = 4.25; % left servo arm length
L2 = 5.75; % left linkage length
L3 = 5.75; % right linkage length
L4 = 4.25; % right servo arm length
L5 = 11/2; % distance from origin to left servo
L6 = 11/2; % distance from origin to right servo

theta_array = [-pi 0.1:0.01:pi-0.1];
phi = 0/180*pi;

for i = 1 : length(theta_array)

theta = theta_array(i);

A(i) = -L3*(-((2*cos(theta)*L1*(sin(phi)*L4-sin(theta)*L1)-2*sin(theta)*L1*(L6 L5-cos(phi)*L4-cos(theta)*L1))/(2*L3*sqrt((L6 L5-cos(phi)*L4-cos(theta)*L1)^2 (sin(phi)*L4-sin(theta)*L1)^2))-((2*sin(theta)*L1*(L6 L5-cos(phi)*L4-cos(theta)*L1)-2*cos(theta)*L1*(sin(phi)*L4-sin(theta)*L1))*(-(L6 L5-cos(phi)*L4-cos(theta)*L1)^2-(sin(phi)*L4-sin(theta)*L1)^2-L3^2 L2^2))/(4*L3*((L6 L5-cos(phi)*L4-cos(theta)*L1)^2 (sin(phi)*L4-sin(theta)*L1)^2)^(3/2)))/sqrt(1-(-(L6 L5-cos(phi)*L4-cos(theta)*L1)^2-(sin(phi)*L4-sin(theta)*L1)^2-L3^2 L2^2)^2/(4*L3^2*((L6 L5-cos(phi)*L4-cos(theta)*L1)^2 (sin(phi)*L4-sin(theta)*L1)^2)))-((cos(theta)*L1)/sqrt((L6 L5-cos(phi)*L4-cos(theta)*L1)^2 (sin(phi)*L4-sin(theta)*L1)^2)-((sin(theta)*L1-sin(phi)*L4)*(2*sin(theta)*L1*(L6 L5-cos(phi)*L4-cos(theta)*L1)-2*cos(theta)*L1*(sin(phi)*L4-sin(theta)*L1)))/(2*((L6 L5-cos(phi)*L4-cos(theta)*L1)^2 (sin(phi)*L4-sin(theta)*L1)^2)^(3/2)))/sqrt(1-(sin(theta)*L1-sin(phi)*L4)^2/((L6 L5-cos(phi)*L4-cos(theta)*L1)^2 (sin(phi)*L4-sin(theta)*L1)^2)))*sin(acos((-(L6 L5-cos(phi)*L4-cos(theta)*L1)^2-(sin(phi)*L4-sin(theta)*L1)^2-L3^2 L2^2)/(2*L3*sqrt((L6 L5-cos(phi)*L4-cos(theta)*L1)^2 (sin(phi)*L4-sin(theta)*L1)^2)))-asin((sin(theta)*L1-sin(phi)*L4)/sqrt((L6 L5-cos(phi)*L4-cos(theta)*L1)^2 (sin(phi)*L4-sin(theta)*L1)^2)));

end

plot(theta_array,A)
  

Результирующий график октавы выглядит следующим образом:

Результат Octave

То же уравнение было скопировано и вставлено из octave в python с заменой ‘^’ на ‘**’, ‘acos’ заменен на ‘arccos’, а ‘asin’ заменен на ‘arcsin’. Был нанесен тот же диапазон тета, при этом phi был равен нулю:

 from pylab import *

# physical setup
L1 = 4.25; # left servo arm length
L2 = 5.75; # left linkage length
L3 = 5.75; # right linkage length
L4 = 4.25; # right servo arm length
L5 = 11.0/2.0; # distance from origin to left servo
L6 = 11.0/2.0; # distance from origin to right servo

theta = arange(-pi 0.1,pi-0.1,0.01);
phi = 0/180.0*pi

def func(theta,phi):

A = -L3*(-((2*cos(theta)*L1*(sin(phi)*L4-sin(theta)*L1)-2*sin(theta)*L1*(L6 L5-cos(phi)*L4-cos(theta)*L1))/(2*L3*sqrt((L6 L5-cos(phi)*L4-cos(theta)*L1)**2 (sin(phi)*L4-sin(theta)*L1)**2))-((2*sin(theta)*L1*(L6 L5-cos(phi)*L4-cos(theta)*L1)-2*cos(theta)*L1*(sin(phi)*L4-sin(theta)*L1))*(-(L6 L5-cos(phi)*L4-cos(theta)*L1)**2-(sin(phi)*L4-sin(theta)*L1)**2-L3**2 L2**2))/(4*L3*((L6 L5-cos(phi)*L4-cos(theta)*L1)**2 (sin(phi)*L4-sin(theta)*L1)**2)**(3/2)))/sqrt(1-(-(L6 L5-cos(phi)*L4-cos(theta)*L1)**2-(sin(phi)*L4-sin(theta)*L1)**2-L3**2 L2**2)**2/(4*L3**2*((L6 L5-cos(phi)*L4-cos(theta)*L1)**2 (sin(phi)*L4-sin(theta)*L1)**2)))-((cos(theta)*L1)/sqrt((L6 L5-cos(phi)*L4-cos(theta)*L1)**2 (sin((phi)*L4-sin(theta)*L1)**2)-((sin(theta)*L1-sin(phi)*L4)*(2*sin(theta)*L1*(L6 L5-cos(phi)*L4-cos(theta)*L1)-2*cos(theta)*L1*(sin(phi)*L4-sin(theta)*L1)))/(2*((L6 L5-cos(phi)*L4-cos(theta)*L1)**2 (sin(phi)*L4-sin(theta)*L1)**2)**(3/2)))/sqrt(1-(sin(theta)*L1-sin(phi)*L4)**2/((L6 L5-cos(phi)*L4-cos(theta)*L1)**2 (sin(phi)*L4-sin(theta)*L1)**2)))*sin(arccos((-(L6 L5-cos(phi)*L4-cos(theta)*L1)**2-(sin(phi)*L4-sin(theta)*L1)**2-L3**2 L2**2)/(2*L3*sqrt((L6 L5-cos(phi)*L4-cos(theta)*L1)**2 (sin(phi)*L4-sin(theta)*L1)**2)))-arcsin((sin(theta)*L1-sin(phi)*L4)/sqrt((L6 L5-cos(phi)*L4-cos(theta)*L1)**2 (sin(phi)*L4-sin(theta)*L1)**2)))

return A

f = figure();
a = f.add_subplot(111);

a.plot(theta,func(theta,phi))

ginput(1, timeout=-1); # wait for user to click so we dont lose the plot
  

Результат Python выглядит следующим образом:
Результат Python

Я не могу определить, что вызывает различия, есть идеи?

Комментарии:

1. Эти функции являются упрощенными версиями исходной функции? Вау. Есть ли шанс, что вы могли бы отрезать одинаковые куски от обоих кусков по одному и попытаться найти что-нибудь еще поменьше ? 🙂

2. Учитывая сложность функции, может ли это быть проблемой различной точности с плавающей запятой и / или ошибок округления? Вы пробовали отображать меньшие части функции, чтобы сузить причину?

3. Это упрощено в том смысле, что весь посторонний код был удален, чтобы упростить проблему для гуру переполнения стека.

4. Что касается ошибок с плавающей запятой / округления, я полагаю , что octave и python используют одинаковый размер (двойная точность) при выполнении вычислений. Честно говоря, это сложное уравнение, и я не отвергаю ошибку плавающей точности / округления / порядка операций

5. tl; dr в python 3/2 = 1, в matlab 3/2 = 1,5

Ответ №1:

Попробуйте from __future__ import division устранить ошибки, возникающие при разделении по этажам.

Комментарии:

1. Ура! Спасибо! похоже, это исправлено. Есть ли какие-либо другие математические ошибки, на которые я должен обратить внимание?

2. @Inverse_Jacobian: Если этот ответ решает вашу проблему, вы должны принять его (нажмите на галочку рядом с ним).