Why is my model not converging with ode45 solver ?
I am trying to simulate an electronical device that can be modeled by a mass-spring-damper system with an additional non-linear force. The equation at the equilibrium for the system is the following :
The goal of my MATLAB code is to solve this equation for $x$ and find the equilibrium. For this purpose, I use the `ode45` function like this (all coefficient are defined in my code but not shown here) :
x0 = 0;
[t, x] = ode45(@(t, x) odefun(t, x, eps0, a, V0, B, d0, K), tspan, x0);
function dxdt = odefun(t, x, eps, a, V, B, d0, K)
dxdt = ((eps * a * V^2) ./ (2 * B * (d0 – x).^2)) + ((K / B) .* x);
end
The equation is good according to my teachers and several papers but the solver never converges and I can’t see why. All coefficient are defined according to the dimensions of a capacitive micromachined ultrasound transducer. The solution of ode45 is the following :
This is obviously wrong because the dimensions of a cmut are bellow millimeter.
Can you see any mistakes in the way I use the ode45 solver ?
Full code :
% DIMENSIONS DE LA MEMBRANE
e = 500e-9; % epaisseur [m]
r = 20e-6; % rayon [m]
d0 = 550e-9; % epaisseur de cavite [m]
a = pi * r^2; % surface de la membrane [m^2]
% PARAMETRES MECANIQUES
E = 200e9; % module d’Young du SiN [Pa]
nu = 0.25; % coefficient de poisson du SiN
eta = 18.5e-6; % viscosité dynamique de l’air [Pa.s]
p0 = 1e5; % pression exterieure [Pa]
rhoSiN = 3170; % densite du SiN [Kg/m^3]
m = rhoSiN * a * e; % masse membrane [Kg]
K = (16 * E * e^3) / (3 * (1 – nu^2) * r^2); % raideur [N/m]
B = (eta * pi * r^2) / e; % amortissement [N.s/m]
% PARAMETRES ELECTRIQUES
eps0 = 8.85e-12; % permittivité diélectrique du vide [F/m]
V0 = 10; % tension de polarisation [V]
% Conditions initiales et plage de temps
x0 = 0;
ti = 0; tf = 1e-3;
dt = 1e-6;
tspan = linspace(ti, tf, 1/dt);
% Résolution par RK4
[t, x] = ode45(@(t, x) odefun(t, x, eps0, a, V0, B, d0, K), tspan, x0);
plot(t,x,’-‘)
function dxdt = odefun(t, x, eps, a, V, B, d0, K)
dxdt = ((eps * a * V^2) ./ (2 * B * (d0 – x).^2)) + ((K / B) .* x);
end
References :
– Y. Wang, L. -M. He, Z. Li, W. Xu and J. Ren, "A computationally efficient nonlinear dynamic model for cMUT based on COMSOL and MATLAB/Simulink"
– T. Merrien, A. Boulmé and D. Certon, "Lumped-Parameter Equivalent Circuit Modeling of CMUT Array Elements"
I tried several solver from MATLAB and implemented my own Runge-Kutta algorithm based on the Wikipedia example. I also verified the coefficients according to my references.I am trying to simulate an electronical device that can be modeled by a mass-spring-damper system with an additional non-linear force. The equation at the equilibrium for the system is the following :
The goal of my MATLAB code is to solve this equation for $x$ and find the equilibrium. For this purpose, I use the `ode45` function like this (all coefficient are defined in my code but not shown here) :
x0 = 0;
[t, x] = ode45(@(t, x) odefun(t, x, eps0, a, V0, B, d0, K), tspan, x0);
function dxdt = odefun(t, x, eps, a, V, B, d0, K)
dxdt = ((eps * a * V^2) ./ (2 * B * (d0 – x).^2)) + ((K / B) .* x);
end
The equation is good according to my teachers and several papers but the solver never converges and I can’t see why. All coefficient are defined according to the dimensions of a capacitive micromachined ultrasound transducer. The solution of ode45 is the following :
This is obviously wrong because the dimensions of a cmut are bellow millimeter.
Can you see any mistakes in the way I use the ode45 solver ?
Full code :
% DIMENSIONS DE LA MEMBRANE
e = 500e-9; % epaisseur [m]
r = 20e-6; % rayon [m]
d0 = 550e-9; % epaisseur de cavite [m]
a = pi * r^2; % surface de la membrane [m^2]
% PARAMETRES MECANIQUES
E = 200e9; % module d’Young du SiN [Pa]
nu = 0.25; % coefficient de poisson du SiN
eta = 18.5e-6; % viscosité dynamique de l’air [Pa.s]
p0 = 1e5; % pression exterieure [Pa]
rhoSiN = 3170; % densite du SiN [Kg/m^3]
m = rhoSiN * a * e; % masse membrane [Kg]
K = (16 * E * e^3) / (3 * (1 – nu^2) * r^2); % raideur [N/m]
B = (eta * pi * r^2) / e; % amortissement [N.s/m]
% PARAMETRES ELECTRIQUES
eps0 = 8.85e-12; % permittivité diélectrique du vide [F/m]
V0 = 10; % tension de polarisation [V]
% Conditions initiales et plage de temps
x0 = 0;
ti = 0; tf = 1e-3;
dt = 1e-6;
tspan = linspace(ti, tf, 1/dt);
% Résolution par RK4
[t, x] = ode45(@(t, x) odefun(t, x, eps0, a, V0, B, d0, K), tspan, x0);
plot(t,x,’-‘)
function dxdt = odefun(t, x, eps, a, V, B, d0, K)
dxdt = ((eps * a * V^2) ./ (2 * B * (d0 – x).^2)) + ((K / B) .* x);
end
References :
– Y. Wang, L. -M. He, Z. Li, W. Xu and J. Ren, "A computationally efficient nonlinear dynamic model for cMUT based on COMSOL and MATLAB/Simulink"
– T. Merrien, A. Boulmé and D. Certon, "Lumped-Parameter Equivalent Circuit Modeling of CMUT Array Elements"
I tried several solver from MATLAB and implemented my own Runge-Kutta algorithm based on the Wikipedia example. I also verified the coefficients according to my references. I am trying to simulate an electronical device that can be modeled by a mass-spring-damper system with an additional non-linear force. The equation at the equilibrium for the system is the following :
The goal of my MATLAB code is to solve this equation for $x$ and find the equilibrium. For this purpose, I use the `ode45` function like this (all coefficient are defined in my code but not shown here) :
x0 = 0;
[t, x] = ode45(@(t, x) odefun(t, x, eps0, a, V0, B, d0, K), tspan, x0);
function dxdt = odefun(t, x, eps, a, V, B, d0, K)
dxdt = ((eps * a * V^2) ./ (2 * B * (d0 – x).^2)) + ((K / B) .* x);
end
The equation is good according to my teachers and several papers but the solver never converges and I can’t see why. All coefficient are defined according to the dimensions of a capacitive micromachined ultrasound transducer. The solution of ode45 is the following :
This is obviously wrong because the dimensions of a cmut are bellow millimeter.
Can you see any mistakes in the way I use the ode45 solver ?
Full code :
% DIMENSIONS DE LA MEMBRANE
e = 500e-9; % epaisseur [m]
r = 20e-6; % rayon [m]
d0 = 550e-9; % epaisseur de cavite [m]
a = pi * r^2; % surface de la membrane [m^2]
% PARAMETRES MECANIQUES
E = 200e9; % module d’Young du SiN [Pa]
nu = 0.25; % coefficient de poisson du SiN
eta = 18.5e-6; % viscosité dynamique de l’air [Pa.s]
p0 = 1e5; % pression exterieure [Pa]
rhoSiN = 3170; % densite du SiN [Kg/m^3]
m = rhoSiN * a * e; % masse membrane [Kg]
K = (16 * E * e^3) / (3 * (1 – nu^2) * r^2); % raideur [N/m]
B = (eta * pi * r^2) / e; % amortissement [N.s/m]
% PARAMETRES ELECTRIQUES
eps0 = 8.85e-12; % permittivité diélectrique du vide [F/m]
V0 = 10; % tension de polarisation [V]
% Conditions initiales et plage de temps
x0 = 0;
ti = 0; tf = 1e-3;
dt = 1e-6;
tspan = linspace(ti, tf, 1/dt);
% Résolution par RK4
[t, x] = ode45(@(t, x) odefun(t, x, eps0, a, V0, B, d0, K), tspan, x0);
plot(t,x,’-‘)
function dxdt = odefun(t, x, eps, a, V, B, d0, K)
dxdt = ((eps * a * V^2) ./ (2 * B * (d0 – x).^2)) + ((K / B) .* x);
end
References :
– Y. Wang, L. -M. He, Z. Li, W. Xu and J. Ren, "A computationally efficient nonlinear dynamic model for cMUT based on COMSOL and MATLAB/Simulink"
– T. Merrien, A. Boulmé and D. Certon, "Lumped-Parameter Equivalent Circuit Modeling of CMUT Array Elements"
I tried several solver from MATLAB and implemented my own Runge-Kutta algorithm based on the Wikipedia example. I also verified the coefficients according to my references. differential equations, physics, matlab, model, ode45, ode MATLAB Answers — New Questions
