How can I divide by the function variable in a Simulink equation?

I want to solve this fairly simple second order differential equation using Simulink:

Which I rewrote as:

The problem is that is another function of r, and because of that I cannot divide by a constant as is normally done when solving differential equations in Simulink.

How can I divide by r in my Simulink model?

[Pat Gipper]

The output of the second integrator should be the solution to your second order equation. Why not feed it into the divide? One word of warning is that Simulink will be solving with a starting value of 0, so you will be dividing by 0 right from the start. To avoid this error try feeding a non-zero initial condition into the second integrator.

[Louis]

Wouldn't feeding it into the divide mean dividing by $v$ then?

Answer 1

[Sam Chak]

Hi @Louis

In the 2nd-order ODE, the independent variable 'r' can represent the radius of the system or another relevant parameter. However, in Simulink, the simulation runs based on the flow of time 't', which is the fundamental independent variable in many dynamic systems. To account for this, you can configure 1/r as if it were 1/t.

This approach, however, may introduce a division-by-zero issue if t starts from 0. To address the potential singularity, you can utilize the Ramp block and configure it accordingly. By setting the Slope to 1, the independent variable r will propagate linearly with respect to Simulink time. Additionally, selecting a positive Initial Output (which corresponds to r) will help avoid the singularity. In the Model Configuration Parameters menu, you can set the Start time to be equal to the Initial Output value.

To verify the results obtained in Simulink, a simple MATLAB code. Comparing the outputs from both the Simulink model and the MATLAB code will help ensure the consistency of your implementation.

k       = 1;
mu      = 1;
ode     = @(r, v) [v(2);
                   k/mu - (1/r)*v(2)];
rspan   = [1, 6];   % r starts from 1 and ends at 6
v0      = [1; 0];   % initial values: v(1) = 1, v'(1) = 0
[r, v]  = ode45(ode, rspan, v0);
plot(r, v(:,1)), grid on, xlabel r, ylabel v, ylim([0 10])

[Louis]

Thanks a lot, this works.