필터 지우기
필터 지우기

Convert symbolic system of differential Equations to numerical for ODE45

조회 수: 2 (최근 30일)
Marius
Marius 2016년 3월 8일
댓글: Star Strider 2016년 3월 8일
Hi all. I have to solve a FEM Problem. For this I created a function which computes the needed matrices and so the set of differential equations. The output of this function is a symbolc vector dy which contains the set of differential equations which are functions of y already in state-space representation. It looks like:
dy = [y5*cos(y3)....;y3*sin(y3)...;...]
I want to read this symbolic vector into a function file say "cylinder_DGL". I already made the computation for an easy problem where the set of differential equations is computed in advance with Mathematica. It looks like:
function dy = cylinder_DGL(t,y,time_r,dt,ug,vg,Sys)
% CYLINDER_DGL evaluates the equation of motion for the rocking and rolling
% cylinder for given time instants t, excitation ug and vg and system
% properties R and H
% Save geometry properties into new variables
r = Sys.r;
h = Sys.h;
g=9.81;
% Interpolate the excitation ug and vg for the time instant t
i = floor(t/dt);
if i < length(time_r)
u = ug(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(ug(i+1)-ug(i));
v = vg(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(vg(i+1)-vg(i));
else
u = 0;
v = 0;
end
% Evaluate the equation of motion with help of the state space
% formulation
% State space formulation: row 1
dy(1,1)= y(2);
% State space formulation: row 2
dy(2,1)= ((4/3*h^2 - 5/4*r^2) * y(4)^2 * cos(y(1)) * sin(y(1))...
+ 3/2*r^2 * y(4)^2 * sin(y(1)) + h*r * y(4)^2 * (1+cos(y(1)))...
- 2*h*r * y(4)^2 * cos(y(1))^2 - r*g * cos(y(1)) + h*g * sin(y(1))...
- h*u * cos(y(3))*cos(y(1)) - r*u * cos(y(3))*sin(y(1))...
- h*v * cos(y(1))*sin(y(3)) - r*v * sin(y(3))*sin(y(1)))...
/(5/4*r^2 + 4/3*h^2);
% State space formulation: row 3
dy(3,1)= y(4);
% State space formulation: row4
dy(4,1)= (-3*r^2 * y(2)*y(4)*sin(y(1))...
+ (5/2*r^2-8/3*h^2)*y(2)*y(4)*sin(y(1))*cos(y(1))...
- 4*h*r*y(2)*y(4)*(sin(y(1))^2+cos(y(1))-1)...
- r*v*cos(y(3))*(1-cos(y(1))) - r*u*sin(y(3))*(1-cos(y(1)))...
- h*v*sin(y(1))*cos(y(3)) + h*u*sin(y(3))*sin(y(1)))...
/ ((4/3*h^2 - 5/4*r^2)*sin(y(1))^2 + 3*r^2*(1-cos(y(1)))...
+ 2*h*r*sin(y(1))*(1-cos(y(1))));
end
This can then be started from the script with:
refine = 4;
options = odeset('Events',@events_overturn,'Refine', refine,...
'RelTol', 1e-9,'absTol',1e-9*[1 1 1 1]);
% Solve the equations of motion
[t,y,te,ye,ie] = ode45(@(t,y) ...
cylinder_DGL(t,y,time_r,dt,ug,vg,Sys),t_3D,IC,options);
How can I read in the symbolic vector dy into my cylinder_DGL file and convert it so that I can use it to be started from the script as for the easy problem?
  댓글 수: 1
Star Strider
Star Strider 2016년 3월 8일
I don’t see that you’re using the Symbolic Math Toolbox anywhere.
In any event, two functions that can help turn a symbolic expression for your ODE into code you can run with the ODE solvers are odeToVectorField and matlabFunction.

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

답변 (0개)

카테고리

Help CenterFile Exchange에서 Symbolic Math Toolbox에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by