2011-11-07 8 views
8

私はオクターブで作成したコードをpylabに移植しています。移植された方程式の1つは、それがオクターブで行うよりもPythonで劇的に異なる結果をもたらします。同じ方程式、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 result

同じ方程式がオクターブからコピーされ、 '^'が '**'、 'acos'が 'arccos'、 'asin'が 'arcsin'に置き換えられて、Pythonにコピーされました。シータの同じ範囲がゼロで開催されたファイでプロットした。

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 result

私は違いの原因を特定カント、任意のアイデア?

+1

これらの関数は、元の関数の_簡略化されたバージョンですか?ワオ。あなたは一度に両方の作品から同じチャンクを取り除くことができ、さらに小さなものをまだ見つけようとするとどんなチャンスですか? :) – sarnold

+0

関数の複雑さを考えると、異なる浮動小数点精度や丸め誤差の問題になる可能性がありますか?原因を絞り込むために関数の小さな部分をプロットしようとしましたか? –

+0

スタックオーバフローGurusの問題を単純化するために、無関係なコードがすべて取り出されているという意味では単純化されています。 –

答えて

12

from __future__ import divisionを試して、床の分割に起因するエラーを排除してください。

+0

Huzzah!ありがとうございました!それはそれを修正したようだ。私が目を向けなければならない他の数学的な問題がありますか? –

+0

@Inverse_Jacobian:この回答があなたの問題を解決したら、それを受け入れるべきです(それによってチェックマークをクリックしてください)。 –

関連する問題