call large equation to ode solver

5 Jun 2018
1 Answer
Answer Accepted
Updated 5 Jun 2018
7 Views (30 days)

0 votes

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

Sign in to answer this question.

 Accepted Answer

Torsten
Torsten on 5 Jun 2018

0 votes

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 Comments
Show 3 older comments Hide 3 older comments

shahin hashemi
shahin hashemi on 5 Jun 2018
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 on 5 Jun 2018
i should mention that GG and G are the same thank you again
Torsten
Torsten on 5 Jun 2018
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 on 5 Jun 2018
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 on 5 Jun 2018
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.

Sign in to comment.

More Answers (0)

Categories

Find more on Numeric Solvers in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!