% Daniel Kawano, Rose-Hulman Institute of Technology % Last modified: Feb 23, 2016 clear all close all clc % Specify symbolic parameters and variables that are functions of the % directed curve's convected coordinate xi: syms psi(xi) theta(xi) phi(xi) syms nu1(xi) nu2(xi) nu3(xi) % (1) Relate the reference triad {D1,D2,D3} = {E1,E2,E3} to the triad % {d1,d2,d3} using a 3-2-3 set of Euler angles: P1 = [cos(psi), sin(psi), 0; -sin(psi), cos(psi), 0; 0, 0, 1]; P2 = [cos(theta), 0, -sin(theta); 0, 1, 0; sin(theta), 0, cos(theta)]; P3 = [cos(phi), sin(phi), 0; -sin(phi), cos(phi), 0; 0, 0, 1]; P = simplify(P3*P2*P1); % (2) Compute the axial vector of strains, nu, which yields a set of % ODEs governing the evolution of the Euler angles for the triad % {d1,d2,d3} along the directed curve: Ptensor = transpose(P); Ktensor = simplify(transpose(Ptensor)*diff(Ptensor, xi)); Ktensor32 = [0, 0, 1]*Ktensor*[0, 1, 0]'; Ktensor13 = [1, 0, 0]*Ktensor*[0, 0, 1]'; Ktensor21 = [0, 1, 0]*Ktensor*[1, 0, 0]'; nu = [Ktensor32; Ktensor13; Ktensor21]; % (3) Manipulate the system of ODEs into a form suitable for numerical % integration: % Compile the state equations and arrange the state variables: StateEqns = [nu1; nu2; nu3] == nu; StateVars = [psi; theta; phi]; % Express the state equations in mass-matrix form, M(xi,Y)*Y'(xi) = % F(xi,Y): [Msym, Fsym] = massMatrixForm(StateEqns, StateVars); Msym = simplify(Msym) Fsym = simplify(Fsym) % Convert M(xi,Y) and F(xi,Y) to symbolic function handles with the input % parameters specified: M = odeFunction(Msym, StateVars, nu1(xi), nu2(xi), nu3(xi)); F = odeFunction(Fsym, StateVars, nu1(xi), nu2(xi), nu3(xi)); % Save M(xi,Y) and F(xi,Y): save directed_curve_ODEs.mat Msym Fsym StateVars M F