Convert symbolic system of differential Equations to numerical for ODE45
1 回表示 (過去 30 日間)
古いコメントを表示
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
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 Center および File Exchange で Symbolic Math Toolbox についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!