Why am i getting the error 'The first input to exist must be a string scalar or character vector' while solving 4 ode using ode45 solver ??
10 views (last 30 days)
Show older comments
I have four differential equation. In these four equation u, s, p, r is the function of t. I solve these equation using ode45 solver but i am getting the error massage. Code is: y0 is the initial condition and its [u0,s0,p0,r0]. Please give me some solution for this error. Thank you for help
clc
clear
close all
b = deg2rad(45);
tspan = [0 pi];
y0 = [0.10 3.14e-06 1.4706 145.60];
syms u(t) s(t) p(t) r(t) t
eq01 = u*diff(s,t) + s*diff(u,t) == (sin(b)^3)/(4*((cos(b)*cos(t)).^2));
eq02 = (u.^2)*diff(s,t)+(2*s*u)*diff(u,t) == cos(b)*(sin(b)^3)/(4*((cos(b)*cos(t)).^2));
eq03 = (u.^2*s)*(1+diff(p,t))*sin(p) == (2*r/52) - (sin(p).^2)*(sin(b)^3)/(4*((cos(b)*cos(t)).^2));
eq04 = diff(r,t) == r*cot(p);
eqs = [eq01, eq02, eq03, eq04];
m = ode45(eqs,tspan,y0);
0 Comments
Answers (2)
Torsten
on 17 Aug 2023
Edited: Torsten
on 17 Aug 2023
ode45 is a numerical integrator. It doesn't accept symbolic inputs.
b = deg2rad(45);
tspan = [0 pi];
y0 = [0.10 3.14e-06 1.4706 145.60];
[T,Y] = ode15s(@(t,y)fun(t,y,b),tspan,y0);
plot(T,Y(:,1))
function dydt = fun(t,y,b)
s = y(1);
u = y(2);
p = y(3);
r = y(4);
M = [u s 0 0;u^2 2*s*u 0 0;0 0 u^2*s*sin(p) 0;0 0 0 1];
rhs = [(sin(b)^3)/(4*((cos(b)*cos(t)).^2));cos(b)*(sin(b)^3)/(4*((cos(b)*cos(t)).^2));(2*r/52) - (sin(p).^2)*(sin(b)^3)/(4*((cos(b)*cos(t)).^2))-u^2*s*sin(p);r*cot(p)];
dydt = M\rhs;
end
2 Comments
Steven Lord
on 17 Aug 2023
b = deg2rad(45);
tspan = [0 pi];
y0 = [0.10 3.14e-06 1.4706 145.60];
[T,Y] = ode15s(@(t,y)fun(t,y,b),tspan,y0);
plot(T,Y(:,1))
% Check reciprocal condition number of M at initial conditions
y = [0.10 3.14e-06 1.4706 145.60];
s = y(1);
u = y(2);
p = y(3);
r = y(4);
format longg
M = [u s 0 0;u^2 2*s*u 0 0;0 0 u^2*s*sin(p) 0;0 0 0 1]
rcond(M)
function dydt = fun(t,y,b)
s = y(1);
u = y(2);
p = y(3);
r = y(4);
M = [u s 0 0;u^2 2*s*u 0 0;0 0 u^2*s*sin(p) 0;0 0 0 1];
rhs = [(sin(b)^3)/(4*((cos(b)*cos(t)).^2));cos(b)*(sin(b)^3)/(4*((cos(b)*cos(t)).^2));(2*r/52) - (sin(p).^2)*(sin(b)^3)/(4*((cos(b)*cos(t)).^2))-u^2*s*sin(p);r*cot(p)];
dydt = M\rhs;
end
Seems to work okay to me (modulo the repeated warnings that you have singular or near-singular matrices involved in your rhs.) Looks like your M matrix, constructed using the values of your initial conditions, is pretty badly conditioned (and I'm guessing element (3, 3) only gets smaller as the ODE solver progresses.)
If you're seeing different behavior please show us the exact code you're running and the full and exact text of any warning and/or error messages you receive (all the text displayed in orange and/or red in the Command Window.)
Though since b only appears in fun as the input to sin or cos I'd skip the deg2rad call (converting 45 degrees into radians) and use the degree-based trig functions sind and cosd instead of the radian-based sin and cos.
Torsten
on 17 Aug 2023
Can you resolve this issue??
No. I just implemented your equations correctly, but you are the only person who knows what they model and thus where the problem integrating them may stem from.
Walter Roberson
on 17 Aug 2023
Here with mechanically generated numeric functions based on your symbolic code. If this does not work then either
- your symbolic expressions are wrong; or
- your initial conditions are wrong; or
- your system is too unstable for numeric computations
b = deg2rad(45);
tspan = [0 pi];
y0 = [0.10 3.14e-06 1.4706 145.60];
syms u(t) s(t) p(t) r(t) t
eq01 = u*diff(s,t) + s*diff(u,t) == (sin(b)^3)/(4*((cos(b)*cos(t)).^2));
eq02 = (u.^2)*diff(s,t)+(2*s*u)*diff(u,t) == cos(b)*(sin(b)^3)/(4*((cos(b)*cos(t)).^2));
eq03 = (u.^2*s)*(1+diff(p,t))*sin(p) == (2*r/52) - (sin(p).^2)*(sin(b)^3)/(4*((cos(b)*cos(t)).^2));
eq04 = diff(r,t) == r*cot(p);
eqs = [eq01, eq02, eq03, eq04];
vars = [u(t), s(t), p(t), r(t)];
[eqs,vars] = reduceDifferentialOrder(eqs, vars)
[M,F] = massMatrixForm(eqs,vars)
f = M\F
odefun = odeFunction(f,vars)
[T,Y] = ode15s(odefun, tspan, y0);
plot(T, Y)
3 Comments
See Also
Categories
Find more on Function Creation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!