0 개 추천

dear all
i want solve equation with using ode113
clc
clear all
tic
initialcond=[-0.001,-0.001,-0.001,-0.001,-0.001,-0.001]
[T X]=ode113(@solver2link,[0 40],initialcond);
toc
this is my ode file and the following file is my equation
function dy = shebhestatici1(t,y)
X1=y(1,1);
X2=y(2,1);X3=y(3,1);X4=y(4,1);X5=y(5,1);X6=y(6,1);
M
G
xdd=inv(M)*(-G);
dy=xdd;
end
i calculate M and G from different file
previously M and G are short equation and i could copy theme directly to my ode solver
but now after some changes it become large file and i cant copy them
so my question is how could i call them to my ode solver for example :
function dy = shebhestatici1(t,y)
X1=y(1,1);
X2=y(2,1);X3=y(3,1);X4=y(4,1);X5=y(5,1);X6=y(6,1);
function M G
xdd=inv(M)*(-G);
dy=xdd;
end
i know its wrong but just to clear my point
thank you all

이 질문에 답변하려면 로그인하십시오.

 채택된 답변

Torsten
Torsten 2018년 6월 5일

0 개 추천

function dy = shebhestatici1(t,y)
X1=y(1,1);
X2=y(2,1);
X3=y(3,1);
X4=y(4,1);
X5=y(5,1);
X6=y(6,1);
M = function_where_you_supply_M(X1,X2,X3,X4,X5,X6);
G = function_where_you_supply_G(X1,X2,X3,X4,X5,X6);
xdd=-M\G;
dy=xdd;
end

댓글 수: 5
이전 댓글 3개 표시 이전 댓글 3개 숨기기

shahin hashemi
shahin hashemi 2018년 6월 5일
ty mr torsten for your answer but i cant understand what do you mean by the where_you_supply_M
this is code where i calculate M and G
clc
clear all
m=7.782e-4;
jx=2.011E-14;
E=2.1e11;
O=0.008;
jz=1.149e-13;
G=8e10;
L=0.09;
g=9.81;
N=4;
T1=0;
T2=0;
T3=0;
Ixx=9.821e-8;
Izz=1.654e-7;
% syms L m jx E O jz G L g Ixx Izz T1 T2 T3
% syms T1 T2 T3
xd = sym('xd', [1 3*N], 'real');
xdd = sym('xdd', [1 3*N], 'real');
x = sym(zeros(1,3*N));%preallocation
syms t
for i=1:3*N
x(i)=sym(sprintf('x%d(t)',i), 'real');
end
rl(:,:,1)=[1*O;0;0;];
rl(:,:,2)=[(-1/2)*O;(sqrt(3))*O/2;0];
rl(:,:,3)=[(-1/2)*O;-(sqrt(3))*O/2;0];
for j=2:N+1
for i=j-1
f1(:,:,j)=(1-cos(x(3*i-2)))/x(3*i-2);
g1(:,:,j)=-(1-cos(x(3*i-1)))/x(3*i-1);
f2(:,:,j)=sin(x(3*i-2))/x(3*i-2);
g2(:,:,j)=sin(x(3*i-1))/x(3*i-1);
f1d(:,:,j)=(sin(x(3*i-2))*x(3*i-2)+cos(x(3*i-2))-1)/(x(3*i-2)^2);
g1d(:,:,j)=-(sin(x(3*i-1))*x(3*i-1)+cos(x(3*i-1))-1)/(x(3*i-1)^2);
f2d(:,:,j)=(cos(x(3*i-2))*x(3*i-2)-sin(x(3*i-2)))/(x(3*i-2)^2);
g2d(:,:,j)=(cos(x(3*i-1))*x(3*i-1)-sin(x(3*i-1)))/(x(3*i-1)^2);
pl(:,:,j)=L*[f1(:,:,j) f2(:,:,j)*g1(:,:,j) f2(:,:,j)*g2(:,:,j)]';
Rl(:,:,j)=[cos(x(3*i))*cos(x(3*i-2)) -sin(x(3*i))*cos(x(3*i-2)) sin(x(3*i-2));cos(x(3*i))*sin(x(3*i-2))*sin(x(3*i-1))+sin(x(3*i))*cos(x(3*i-1)) -sin(x(3*i))*sin(x(3*i-2))*sin(x(3*i-1))+cos(x(3*i))*cos(x(3*i-1)) -cos(x(3*i-2))*sin(x(3*i-1));-cos(x(3*i))*sin(x(3*i-2))*cos(x(3*i-1))+sin(x(3*i))*sin(x(3*i-1)) sin(x(3*i))*sin(x(3*i-2))*cos(x(3*i-1))+cos(x(3*i))*sin(x(3*i-1)) cos(x(3*i-2))*cos(x(3*i-1))];
p(:,:,2)=pl(:,:,2);
end
end
R(:,:,2)=Rl(:,:,2);
R(:,:,1)=[1 0 0;0 1 0;0 0 1];
for i=3:N+1
R(:,:,i)=R(:,:,i-1)*Rl(:,:,i);
p(:,:,i)=p(:,:,i-1)+R(:,:,i-1)*pl(:,:,i);
end
for j=2:N+1
for i=j-1
wl(:,:,j)=[sin(x(3*i-2))*(diff(x(3*i), t))+diff(x(3*i-1),t);-cos(x(3*i-2))*sin(x(3*i-1))*(diff(x(3*i), t))+cos(x(3*i-1))*(diff(x(3*i-2), t));cos(x(3*i-2))*cos(x(3*i-1))*(diff(x(3*i), t))+sin(x(3*i-1))*(diff(x(3*i-2), t))];
w(:,:,2)=wl(:,:,2);
end
end
for i=3:N+1
w(:,:,i)= w(:,:,i-1)+R(:,:,i-1)*wl(:,:,i);
end
for j=2:N+1
for i=j-1
pd(:,:,j)=L*[f1d(:,:,j) 0 0;f2d(:,:,j)*g1(:,:,j) f2(:,:,j)*g1d(:,:,j) 0;f2d(:,:,j)*g2(:,:,j) f2(:,:,j)*g2d(:,:,j) 0]*[(diff(x(3*i-2), t));(diff(x(3*i-1), t));(diff(x(3*i), t))];
v(:,:,2)=pd(:,:,2);
end
end
for i=3:N+1
v(:,:,i)=v(:,:,i-1)+cross(w(:,:,i-1),R(:,:,i-1)*pl(:,:,i))+R(:,:,i-1)*pd(:,:,i);
end
for j=2:N+1
wdlcl(:,:,j)=diff(wl(:,:,j),t);
alfa(:,:,2)=wdlcl(:,:,2);
end
for i=3:N+1
alfa(:,:,i)= alfa(:,:,i-1)+cross(w(:,:,i-1),R(:,:,i-1)*wl(:,:,i))+R(:,:,i-1)*wdlcl(:,:,i);
end
for i=2:N+1
for j=1:N
alfa(:,:,i)=subs(alfa(:,:,i),(diff(x(3*j), t)),xd(3*j));
alfa(:,:,i)=subs(alfa(:,:,i),(diff(x(3*j-1), t)),xd(3*j-1));
alfa(:,:,i)=subs(alfa(:,:,i),(diff(x(3*j-2), t)),xd(3*j-2));
alfa(:,:,i)=subs(alfa(:,:,i),(diff(x(3*j), t,t)),xdd(3*j));
alfa(:,:,i)=subs(alfa(:,:,i),(diff(x(3*j-1), t,t)),xdd(3*j-1));
alfa(:,:,i)=subs(alfa(:,:,i),(diff(x(3*j-2), t,t)),xdd(3*j-2));
end
end
for j=2:N+1
pdd(:,:,j)=diff(pd(:,:,j),t);
a(:,:,2)=pdd(:,:,2);
end
for i=3:N+1
a(:,:,i)= a(:,:,i-1)+cross(diff(w(:,:,i-1),t),R(:,:,i-1)*pl(:,:,i))+cross(2*w(:,:,i-1),R(:,:,i-1)*pd(:,:,i))+cross(w(:,:,i-1),cross(w(:,:,i-1),R(:,:,i-1)*pl(:,:,i-1)))+R(:,:,i-1)*pdd(:,:,i);
end
for i=2:N+1
for j=1:N
a(:,:,i)=subs(a(:,:,i),(diff(x(3*j), t)),xd(3*j));
a(:,:,i)=subs(a(:,:,i),(diff(x(3*j-1), t)),xd(3*j-1));
a(:,:,i)=subs(a(:,:,i),(diff(x(3*j-2), t)),xd(3*j-2));
a(:,:,i)=subs(a(:,:,i),(diff(x(3*j), t,t)),xdd(3*j));
a(:,:,i)=subs(a(:,:,i),(diff(x(3*j-1), t,t)),xdd(3*j-1));
a(:,:,i)=subs(a(:,:,i),(diff(x(3*j-2), t,t)),xdd(3*j-2));
w(:,:,i)=subs(w(:,:,i),(diff(x(3*j), t)),xd(3*j));
w(:,:,i)=subs(w(:,:,i),(diff(x(3*j-1), t)),xd(3*j-1));
w(:,:,i)=subs(w(:,:,i),(diff(x(3*j-2), t)),xd(3*j-2));
v(:,:,i)=subs(v(:,:,i),(diff(x(3*j), t)),xd(3*j));
v(:,:,i)=subs(v(:,:,i),(diff(x(3*j-1), t)),xd(3*j-1));
v(:,:,i)=subs(v(:,:,i),(diff(x(3*j-2), t)),xd(3*j-2));
end
end
for i=2:N+1
Fi(:,:,i)=-m*a(:,:,i);
end
for i=2:N+1
Ii(:,:,i)=R(:,:,i)*[Ixx 0 0;0 Ixx 0;0 0 Izz]*transpose(R(:,:,i));
Mi(:,:,i)=-Ii(:,:,i)*alfa(:,:,i)-cross(w(:,:,i),Ii(:,:,i)*w(:,:,i));
end
for i=2:N+1
wk(:,:,i)=jacobian(w(:,:,i),[xd(1),xd(2),xd(3) ,xd(4),xd(5),xd(6),xd(7),xd(8),xd(9),xd(10),xd(11),xd(12)]);% ,xd(13),xd(14),xd(15),xd(16),xd(17),xd(18),xd(19),xd(20),xd(21),xd(22),xd(23),xd(24)]);
end
for i=2:N+1
vk(:,:,i)=jacobian(v(:,:,i),[xd(1),xd(2),xd(3) ,xd(4),xd(5),xd(6),xd(7),xd(8),xd(9),xd(10),xd(11),xd(12)]);%,xd(13),xd(14),xd(15),xd(16),xd(17),xd(18),xd(19),xd(20),xd(21),xd(22),xd(23),xd(24)]);
end
for j=N+1
for i=j-1
Mb(:,:,j)=E*jx*((x(3*i-1))/L)*R(:,:,j-1)*[1;0;0]+E*jx*((x(3*i-2))/L)*R(:,:,j-1)*[0;1;0];
Mt(:,:,j)=-G*jz*(x(3*i))*R(:,3,j)/L;
Me(:,:,j)=Mt(:,:,j)-Mb(:,:,j);
end
end
for j=2:N+1
for i=j-1
Mb(:,:,j)=E*jx*((x(3*i-1))/L)*R(:,:,j-1)*[1;0;0]+E*jx*((x(3*i-2))/L)*R(:,:,j-1)*[0;1;0];
end
end
for j=2:N
for i=j-1
Mt(:,:,j)=-G*jz*(x(3*i))*R(:,3,j)/L;
end
end
for j=2:N
for i=j
Mtt(:,:,j)=G*jz*(x(3*i))*R(:,3,j)/L;
end
end
for j=N+1
Mtt(:,:,j)=[0;0;0];
end
for i=2:N+1
MT(:,:,i)=Mt(:,:,i)+Mtt(:,:,i);
end
for i=2:N
Me(:,:,i)=Mb(:,:,i+1)-Mb(:,:,i)+MT(:,:,i);
end
for i=2:N+1
Fg(:,:,i)=-m*g*[1;0;0];
end
for i=2
ph(:,1,i)=[pl(:,:,i)+R(:,:,i)*rl(:,:,1)-rl(:,:,1)];
ph(:,2,i)=[pl(:,:,i)+R(:,:,i)*rl(:,:,2)-rl(:,:,2)];
ph(:,3,i)=[pl(:,:,i)+R(:,:,i)*rl(:,:,3)-rl(:,:,3)];
end
for i=3:N+1
ph(:,1,i)=[R(:,:,i-1)*pl(:,:,i)-R(:,:,i-1)*rl(:,:,1)+R(:,:,i)*rl(:,:,1)];
ph(:,2,i)=[R(:,:,i-1)*pl(:,:,i)-R(:,:,i-1)*rl(:,:,2)+R(:,:,i)*rl(:,:,2)];
ph(:,3,i)=[R(:,:,i-1)*pl(:,:,i)-R(:,:,i-1)*rl(:,:,3)+R(:,:,i)*rl(:,:,3)];
end
for i=2:N+1
ap(:,:,i)=[sqrt((ph(1,1,i))^2 + (ph(2,1,i))^2 + (ph(3,1,i))^2 );sqrt((ph(1,2,i))^2 + (ph(2,2,i))^2 + (ph(3,2,i))^2);sqrt((ph(1,3,i))^2 + (ph(2,3,i))^2 + (ph(3,3,i))^2)];
f(:,:,i)=[ph(:,1,i)/ap(1,1,i) ph(:,2,i)/ap(2,1,i) ph(:,3,i)/ap(3,1,i) ];
end
for i=N+1
Fc(:,:,i)=[-T1*f(:,1,i) -T2*f(:,2,i) -T3*f(:,3,i) ];
end
for i=2:N
Fc(:,:,i)=[T1*(f(:,1,i+1)-f(:,1,i)) T2*(f(:,2,i+1)-f(:,2,i)) T3*(f(:,3,i+1)-f(:,3,i)) ];
end
for i=2:N+1
fa(:,:,i)=Fc(:,1,i)+Fc(:,2,i)+Fc(:,3,i);
end
for i=2:N+1
Ma(:,:,i)=cross(R(:,:,i)*rl(:,:,1),Fc(:,1,i))+cross(R(:,:,i)*rl(:,:,2),Fc(:,2,i))+cross(R(:,:,i)*rl(:,:,3),Fc(:,3,i));
end
for i=2:N+1
feq(:,:,i)=fa(:,:,i)+Fg(:,:,i)+Fi(:,:,i);
Meq(:,:,i)=Ma(:,:,i)+Me(:,:,i)+Mi(:,:,i);
end
for j=1:(3*N)
for i=2:N+1
FF(j,i-1)= dot(Meq(:,:,i),wk(:,j,i))+dot(feq(:,:,i),vk(:,j,i));
end
end
for j=1:(3*N)
F(j,1)=sum(FF(j,:));
end
for i=1:3*N
M(i,:)=jacobian(F(i,1),[xdd(1),xdd(2),xdd(3) ,xdd(4),xdd(5),xdd(6),xdd(7),xdd(8),xdd(9),xdd(10),xdd(11),xdd(12)]);%,xdd(13),xdd(14),xdd(15),xdd(16),xdd(17),xdd(18),xdd(19),xdd(20),xdd(21),xdd(22),xdd(23),xdd(24)]);
end
% % for i=1:3*N
% % for j=1:N
% % F(i,1)=subs(F(i,1),xd(3*j), 0);
% % F(i,1)=subs(F(i,1),xd(3*j-1), 0);
% % F(i,1)=subs(F(i,1),xd(3*j-2), 0);
% % end
% % end
% % for i=1:3*N
% % M(i,:)=jacobian(F(i,1),[xdd(1),xdd(2),xdd(3) ,xdd(4),xdd(5),xdd(6)]);%,,xdd(7),xdd(8),xdd(9),xdd(10),xdd(11),xdd(12)]);%,xdd(13),xdd(14),xdd(15),xdd(16),xdd(17),xdd(18),xdd(19),xdd(20),xdd(21),xdd(22),xdd(23),xdd(24)]);
% % end
% %
% % test M ba F dar zamane jacobian
GG=F;
for i=1:3*N
for j=1:N
GG(i,1)=subs(GG(i,1),xd(3*j), 0);
GG(i,1)=subs(GG(i,1),xd(3*j-1), 0);
GG(i,1)=subs(GG(i,1),xd(3*j-2), 0);
GG(i,1)=subs(GG(i,1),xdd(3*j), 0);
GG(i,1)=subs(GG(i,1),xdd(3*j-1), 0);
GG(i,1)=subs(GG(i,1),xdd(3*j-2), 0);
end
end
X = sym('X', [1 3*N], 'real');
for i=1:3*N
for j=1:N
GG(i,1)=subs(GG(i,1),x(3*j),X(3*j));
GG(i,1)=subs(GG(i,1),x(3*j-1),X(3*j));
GG(i,1)=subs(GG(i,1),x(3*j-2), X(3*j));
end
end
for i=1:3*N
for j=1:N
M(i,:)=subs(M(i,:),x(3*j),X(3*j));
M(i,:)=subs(M(i,:),x(3*j-1),X(3*j));
M(i,:)=subs(M(i,:),x(3*j-2), X(3*j));
end
end
and i save this file with the name of shebhestatici1 can you tell me that what changes should i do ?
shahin hashemi
shahin hashemi 2018년 6월 5일
i should mention that GG and G are the same thank you again
Torsten
Torsten 2018년 6월 5일
Make your code a function with output argments G and M and call it from where you need M and G.
shahin hashemi
shahin hashemi 2018년 6월 5일
i can change my code to a function but just with one output
function M = shebhestatici1
.
.
.
.
end
how should i do that for 2 ?
maybe like this ?
this doesnt work
function [M G] = shebhestatici1
and for last question should i put smth in front of shebhestatici1 ?
like :
function M = shebhestatici1(x)
Torsten
Torsten 2018년 6월 5일
https://de.mathworks.com/help/matlab/ref/function.html
In your case:
function [M,G] = shebhestatici1
If X1,...,X6 are needed to generate M and G,
function [M,G] = shebhestatici1(X1,X2,X3,X4,X5,X6)
Best wishes
Torsten.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

도움말 센터 및 File Exchange에서 Numeric Solvers에 대해 자세히 알아보기

태그

질문:

shahin hashemi
shahin hashemi
2018년 6월 5일

댓글:

Torsten
Torsten
2018년 6월 5일

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by