ode45 solver argument error

clear all

clc;

syms y(x)

K = 0.1;

A=10;

N=100;

v=diff(y);

Y0 = [y(0)==K, v(0)==0];

tspan = [0 20];

eqn = diff(y,2) == N*(1+v^2)^1.5 + A*(K - y)*(1+v^2)^1.5 - (1+v^2)*v*x^(-1);

V=odeToVectorField(eqn,Y0);

M = matlabFunction(V,'vars',{'x','Y'});

ySol = ode45(M,tspan,Y0);

xValues = linspace(0,20,100);

yValues = deval(ySol,xValues,1);

plot(xValues,yValues)

I had a problem with ySol line in its arguments but I put in everything as i read from matlab informational page but then arithmetical error appeared(probably during tinkering with code i changed something) but I really dont have an idea what causes the main error, its really hard to understand error meaning, I am using odetovectorField for second order differential equation solution, if anyone knows how to successfully solve this diff eq please share. anyway here are errors

Error using ./

Arithmetical expression expected.

Error in symengine>@(x,Y)[Y(2);-(Y(1).*1.0e+1-1.0).*(Y(2).^2+1.0).^(3.0./2.0)+(Y(2).^2+1.0).^(3.0./2.0).*1.0e+2-((Y(2).^2+1.0).*Y(2))./x]

Error in odearguments (line 92)

f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0. this is error i dont understand, is it connected to line 14?

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in (line 14)

ySol = ode45(M,tspan,Y0);

[Dyuman Joshi]

The input y0 in ode45 i.e. the initial conditions are supposed to be numeric, but you have defined them to be symbolic, which is why you get the error.

syms y(x)
K = 0.1;
A = 10;
N = 100;
v = diff(y);
Y0 = [y(0)==K, v(0)==0];
tspan = [0 20];
eqn = diff(y,2) == N*(1+v^2)^1.5 + A*(K - y)*(1+v^2)^1.5 - (1+v^2)*v*x^(-1);
V = odeToVectorField(eqn);
M = matlabFunction(V,'vars',{'x','Y'});

%Supplying numeric values for initial conditions
ySol = ode45(M,tspan,[K; 0]);
xValues = linspace(0,20,100)

xValues = 1×100

0 0.2020 0.4040 0.6061 0.8081 1.0101 1.2121 1.4141 1.6162 1.8182 2.0202 2.2222 2.4242 2.6263 2.8283 3.0303 3.2323 3.4343 3.6364 3.8384 4.0404 4.2424 4.4444 4.6465 4.8485 5.0505 5.2525 5.4545 5.6566 5.8586

yValues = deval(ySol,xValues,1)

yValues = 1×100

0.1000 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

%plot(xValues,yValues)

[Andria]

thank you understood

Answer 1

[Sam Chak]

Hi @Andria

If you start the simulation from , then a singularity will occur due to the division-by-zero term '...*v/x'. However, the system is inherently unstable, as you can see the state v exponentially explodes to a very large value in a very short interval of x.

xspan  = [0.001 0.02];  % <-- start from non-zero x that is relatively close to 0
Y0     = [0.1; 0];
[x, Y] = ode45(@M, xspan, Y0);

plot(x, Y, 'linewidth', 1.5), grid on, 
legend({'$y$', '$v$'}, 'interpreter', 'latex', 'fontsize', 12), legend('boxoff')
xlabel('x'), ylabel('Amplitude')

% the system
function Ydot = M(x, Y)
    % definitions
    Ydot    = zeros(2, 1);
    y       = Y(1);
    v       = Y(2);
    
    % parameters
    K       = 0.1;
    A       = 10;
    N       = 100;
    
    % system
    Ydot(1) =  v;
    Ydot(2) =  N*(1 + v^2)^1.5 + A*(K - y)*(1 + v^2)^1.5 - (1 + v^2)*v/x;
% diff(y,2) == N*(1 + v^2)^1.5 + A*(K - y)*(1 + v^2)^1.5 - (1 + v^2)*v*x^(-1);
end