How to solve a system of coupled first order differential equations using ode45?

1 次查看（过去 30 天）
Md. Golam Zakaria
评论： dpb 2023-10-25
Hello everyone,
I have a system of four coupled first oder differential equation, which needs to be solved. While the equations are long, its pretty straightforward. I have written the following code and getting some error.
function dydt=ccp_ode(t,y)
A_E=.00707;
A_G=0.1650;
e=1.60217663e-19;
n=3.15e16;
l_B=0.1;
epsilon_0=8.85418782e-12;
v_e = 594.0413;
m_e=9.10938e-31;
C_B=2e-9;
m_i=40;
T_e=1.78;
u_B=9.8e3*sqrt(T_e/m_i);
T_ar=300;
p=1;
k_B=1.380649e-23;
n_g=(p/(k_B*T_ar));
I_iP=e*n*u_B*A_E;
I_iG=e*n*u_B*A_G;
L_pl=(l_B*m_e)/(e^2*n*A_E);
K_iz= 2.34e-14*T_e^0.59*exp(-17.44/T_e);
K_ex=2.48e-14*T_e^0.33*exp(-12.78/T_e);
K_el=2.336e-14*T_e^1.609*exp(0.0618*(log(T_e))^2-0.1171*(log(T_e))^3);
Km=K_iz+K_ex+K_el;
v_m=Km*n_g;
v_eff=v_m+(v_e/l_B);
A=2*e*epsilon_0*n*A_E^2;
B=2*e*epsilon_0*n*A_G^2;
C=e*n*v_e*A_E;
D=e*n*v_e*A_G;
f=13.56e6;
omega_1=2*pi*f;
V0=300;
dydt(1)=((-(sqrt(A/y(1)))^-1)*(y(3)+I_iP-(C*exp((-e*y(1))/T_e))));
dydt(2)=((-(sqrt(B/y(2)))^-1)*(-y(3)+I_iG-(D*exp((-e*y(2))/T_e))));
dydt(3)=((L_pl^-1)*((V0*cos(omega_1*t))+y(1)-y(2)+y(4))-(v_eff^y(3)));
dydt(4)=((-1/C_B)*y(3));
end
To call this function I used [t,y]=ode45('ccp_ode', [0 .1], [0 0 2 0]). However I am getting the following error:
Error using odearguments (line 93)
CCP_ODE must return a column vector.
Can anyone please , point me the problem with my code and what is the solution.
Thank You.
  1 个评论
MOSLI KARIM
MOSLI KARIM 2023-10-25
function d
[t,y]=ode45(@ccp_ode, [0 .1], [0 0 2 0])
plot(t,y)
function DYDT=ccp_ode(t,y)
A_E=.00707;
A_G=0.1650;
e=1.60217663e-19;
n=3.15e16;
l_B=0.1;
epsilon_0=8.85418782e-12;
v_e = 594.0413;
m_e=9.10938e-31;
C_B=2e-9;
m_i=40;
T_e=1.78;
u_B=9.8e3*sqrt(T_e/m_i);
T_ar=300;
p=1;
k_B=1.380649e-23;
n_g=(p/(k_B*T_ar));
I_iP=e*n*u_B*A_E;
I_iG=e*n*u_B*A_G;
L_pl=(l_B*m_e)/(e^2*n*A_E);
K_iz= 2.34e-14*T_e^0.59*exp(-17.44/T_e);
K_ex=2.48e-14*T_e^0.33*exp(-12.78/T_e);
K_el=2.336e-14*T_e^1.609*exp(0.0618*(log(T_e))^2-0.1171*(log(T_e))^3);
Km=K_iz+K_ex+K_el;
v_m=Km*n_g;
v_eff=v_m+(v_e/l_B);
A=2*e*epsilon_0*n*A_E^2;
B=2*e*epsilon_0*n*A_G^2;
C=e*n*v_e*A_E;
D=e*n*v_e*A_G;
f=13.56e6;
omega_1=2*pi*f;
V0=300;
dydt(1)=((-(sqrt(A/y(1)))^-1)*(y(3)+I_iP-(C*exp((-e*y(1))/T_e))));
dydt(2)=((-(sqrt(B/y(2)))^-1)*(-y(3)+I_iG-(D*exp((-e*y(2))/T_e))));
dydt(3)=((L_pl^-1)*((V0*cos(omega_1*t))+y(1)-y(2)+y(4))-(v_eff^y(3)));
dydt(4)=((-1/C_B)*y(3));
DYDT=dydt'
end
end

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

dpb
dpb 2023-2-1
...
dydt(1)=((-(sqrt(A/y(1)))^-1)*(y(3)+I_iP-(C*exp((-e*y(1))/T_e))));
dydt(2)=((-(sqrt(B/y(2)))^-1)*(-y(3)+I_iG-(D*exp((-e*y(2))/T_e))));
dydt(3)=((L_pl^-1)*((V0*cos(omega_1*t))+y(1)-y(2)+y(4))-(v_eff^y(3)));
dydt(4)=((-1/C_B)*y(3));
The error message tells you what the problem is -- the above will return a row vector because you didn't preallocate the variable nor expressly define the terms as column vector members nor reshape() the result to ensure the function returns a column vector.
You could have found this easily by just using the debugger to inpsect the result in your function before it returns or looking at what it returned if called from the command line standalone.
To solve it is simple; several ways...
  1. Preallocate the vector before populating it...
dydt=zeros(4,1);
dydt(1)= ....
...
2. Explicitly store into first column...
dydt(1,1)=...
dydt(2,1)=...
...
3. Ensure result is column vector...
...
dydt(1)=...
dydt(2)=...
...
dydt=dydt(:);
return
end
Of these, the cleanest is the first as it avoids the dynamic reallocation that occurs with any of the others without having done so. For only four elements the performance difference will be negligible except the function will be called every solution pass and multiple times for evaluation for every time step so performance here may well be significant. Hence, get into the habit of preallocating arrays where can.
  3 个评论
显示 1更早的评论
MOSLI KARIM
MOSLI KARIM 2023-10-25
编辑：dpb 2023-10-25
my proposed solution
function d
[t,y]=ode45(@ccp_ode, [0 .1], [0 0 2 0])
plot(t,y)
function DYDT=ccp_ode(t,y)
A_E=.00707;
A_G=0.1650;
e=1.60217663e-19;
n=3.15e16;
l_B=0.1;
epsilon_0=8.85418782e-12;
v_e = 594.0413;
m_e=9.10938e-31;
C_B=2e-9;
m_i=40;
T_e=1.78;
u_B=9.8e3*sqrt(T_e/m_i);
T_ar=300;
p=1;
k_B=1.380649e-23;
n_g=(p/(k_B*T_ar));
I_iP=e*n*u_B*A_E;
I_iG=e*n*u_B*A_G;
L_pl=(l_B*m_e)/(e^2*n*A_E);
K_iz= 2.34e-14*T_e^0.59*exp(-17.44/T_e);
K_ex=2.48e-14*T_e^0.33*exp(-12.78/T_e);
K_el=2.336e-14*T_e^1.609*exp(0.0618*(log(T_e))^2-0.1171*(log(T_e))^3);
Km=K_iz+K_ex+K_el;
v_m=Km*n_g;
v_eff=v_m+(v_e/l_B);
A=2*e*epsilon_0*n*A_E^2;
B=2*e*epsilon_0*n*A_G^2;
C=e*n*v_e*A_E;
D=e*n*v_e*A_G;
f=13.56e6;
omega_1=2*pi*f;
V0=300;
dydt(1)=((-(sqrt(A/y(1)))^-1)*(y(3)+I_iP-(C*exp((-e*y(1))/T_e))));
dydt(2)=((-(sqrt(B/y(2)))^-1)*(-y(3)+I_iG-(D*exp((-e*y(2))/T_e))));
dydt(3)=((L_pl^-1)*((V0*cos(omega_1*t))+y(1)-y(2)+y(4))-(v_eff^y(3)));
dydt(4)=((-1/C_B)*y(3));
DYDT=dydt'
end
end
dpb
dpb 2023-10-25
function DYDT=ccp_ode(t,y)
...
dydt(1)=((-(sqrt(A/y(1)))^-1)*(y(3)+I_iP-(C*exp((-e*y(1))/T_e))));
dydt(2)=((-(sqrt(B/y(2)))^-1)*(-y(3)+I_iG-(D*exp((-e*y(2))/T_e))));
dydt(3)=((L_pl^-1)*((V0*cos(omega_1*t))+y(1)-y(2)+y(4))-(v_eff^y(3)));
dydt(4)=((-1/C_B)*y(3));
DYDT=dydt';
end
NOTA BENE: -- The above doesn't preallocate the array and also creates a new temporary variable just to reshape the existing vector. The ' operator also is conjugate transpose; while it should not generate complex values here, it is the .' operator for nonconjugate transpose that is wanted here; it would change the meaning of DYDT from that of dydt if were complex.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

标签

产品


版本

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by