# Using ODE45 to Solve a System of Equations Outputting Linear non-Ideal Results

2 views (last 30 days)
John on 17 Feb 2024
Edited: Torsten on 20 Feb 2024
I am currently modeling a thruster in 1D. And using a set of governing equations found in a paper by Tommaso Misuri, the equations are 3.1 to 3.6. The form I have rearranged the equations in to is:
M(Y) * dY/dz = F(Y)
Where M(Y) is the mass matrix, dY/dz are the collection of diff equations, and F(Y) is the vector of RHS.
The following is the function for solving this system.
function dY_dz = odefcn(z, Y, sigma, rho, R, alpha, muo, j)
v = Y(1);
Tk = Y(2);
B_self = Y(3);
c_v = (3/2)*R; % Specific heat at specific volume (ideal monatomic)
c_p = c_v + R; % Specific heat ratio at specific pressure (ideal monatomic)
% F(Y) Right-Hand Side
rhs_eq = zeros(size(Y));
rhs_eq(1) = (1/rho)*(j*B_self); % (1) Momentum Equation
rhs_eq(2) = (j*v*B_self) + (((j^2)/sigma)*(1 - alpha)); % (2) Energy Equation
rhs_eq(3) = -muo*j; % (3) Maxwell's equation for self induced mag field
% Mass Matrix
M = zeros(length(Y), length(Y));
M(1, 1) = v; % (1) Momentum Equation
M(2, 1) = rho*(v^2); M(2, 2) = rho*v*c_p; % (2) Energy Equation
M(3, 3) = 1; % (3) B_self
dY_dz = M \ rhs_eq;
end
And this is the setup for ODE45:
v_o = velocity;
T_o = Tk;
B_o = B_self;
Y0 = [v_o; T_o; B_o]; % Initial conditions
[z, Y] = ode45(@(z, Y) odefcn(z, Y, sigma, rho, R, alpha, muo, j), ...
z, ... % tspan
Y0); % Initial conditions
And these are the results:
As you can see these are less than ideal, for an electric propulsion thruster velocity increase should be dramatically higher than this (and non-linear). Ideal results should look similar to the paper linked previously.
My initial thoughts as to why they are incorrect is because the function I am using 'reuses' the v, T, and B_self variables, so basically the function does compute a new value but then as soon as it is looped back into the function it just reassigns it back to the initial value. This would explain why the results do not change greatly.
Would this be the reason why my code doesn't work? Or is it a different problem? My knowledge on ODE45 is limited. Any help will be appreciated.
John on 20 Feb 2024
As the paper I'm referencing from states that the continuity equation is: rho*v*A = constant. I just assumed constant is zero. Is this incorrect?
Regarding ode15s and ode23t neither of them worked unfortunately.
Torsten on 20 Feb 2024
Edited: Torsten on 20 Feb 2024
I just assumed constant is zero. Is this incorrect?
rho*v has unit kg/(m^2*s) and is the mass flow per unit area into the domain. Setting it to 0 is a very uninteresting case because your initial value either for rho or v at z=0 had to be 0.
rho(z)*v(z)*A = rho0*v0*A
where rho0, v0 are density and velocity at z = 0.
Thus your first equation must be
rho(z)*v(z) - rho0*v0 = 0
assuming A does not change with z.

Steven Lord on 17 Feb 2024
The form I have rearranged the equations in to is:
M(Y) + dY/dz = F(Y)
Where M(Y) is the mass matrix, dY/dz are the collection of diff equations, and F(Y) is the vector of RHS.
Can you confirm this is the correct form of the equations? When a mass matrix is involved usually it multiplies the derivative, M(Y)*dy/dz = F(Y), rather than adding to it as you've written it. If you are solving the case with a mass matrix multiplying the derivative, rather than trying to handle it yourself in your ODE function I'd use odeset to set the Mass option and pass that options structure into ode45.
If the addition form is correct, your code is incorrect. The right-hand side should be F(Y)-M(Y).
This documentation page works through an example with a mass matrix; since you say your knowledge of ode45 is limited perhaps you could use this as a guide or model for how to solve an ODE with a mass matrix.
John on 19 Feb 2024
Hello Steven,
Thankyou for pointing this out, that was my error, it should've been multiply. I have changed this.
Torsten on 19 Feb 2024
In your code, you multiplied from the very beginning ...

### Categories

Find more on Ordinary Differential Equations in Help Center and File Exchange

R2022b

### Community Treasure Hunt

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

Start Hunting!