solving non-linear ODE

I'm trying to solve ODE using MATLAB(ode45), but not working.

In this case, how can I modify code

Here is existing code

a = theta_m / erf(Z)/(2*sqrt(alpha * t));

c = (1 - theta_m)/(erfc(Z)/(2*sqrt(alpha * t)));

term1 = simplify(a - c);

term2 = (sqrt(pi)*rho*del_H/(2*k*(T1-T0))) * exp((Z)^2/(4*alpha^2*t^2));

dZdt = @(t,Z) term1/term2

tspan = [0 5];

Z0 = 0;

[t,Z] = ode45(dZdt, tspan, Z0);

Error using superiorfloat

Inputs must be floats, namely single or double.

Error in odearguments (line 114)

dataType = superiorfloat(t0,y0,f0);

Error in ode45 (line 107)

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

Error in untitled3 (line 27)

[t,Z] = ode45(dZdt, tspan, Z0);

Thanks for your support

Answer 1

[Sam Chak]

Hi @문기

I didn't modify your equations (except for removing the 'simplify' part), but you need to provide values for the parameters in the 'ode' function. Some division terms encounter a division-by-zero singularity (referred to as Infinity ∞ in layman's terms) when and . Therefore, I adjusted the code by using small positive values for tspan(1) and Z0.

%% Define the differential equation in a function object
function dZdt = ode(t, Z)
    % parameters
    theta_m = 1;
    alpha   = 1;
    rho     = 1;
    del_H   = 1;
    k       = 1;
    T1      = 1;
    T0      = 0;
    
    % terms
    a       = theta_m / erf(Z) / (2*sqrt(alpha*t));
    c       = (1 - theta_m)/(erfc(Z)/(2*sqrt(alpha * t)));
    term1   = a - c;
    term2   = (sqrt(pi)*rho*del_H/(2*k*(T1-T0))) * exp((Z)^2/(4*alpha^2*t^2));
    
    % ODE
    dZdt    = term1/term2;
end

tspan   = [0.0001 5];
Z0      =  0.0001;
[t, Z]  = ode45(@ode, tspan, Z0);
plot(t, Z), grid on, xlabel('t'), ylabel('Z(t)'), title('Time response of Z')

[문기]

the key was the fact that it cannot be divided with 0.

Thank you

Answer 2

[Steven Lord]

This line suggests you're performing operations with symbolic variables.

term1 = simplify(a - c);

If so, you're probably going to need to use matlabFunction or odeFunction to create the anonymous function rather than simply dividing your terms (which won't substitute values into the symbolic expression) or use dsolve.

See the section on solving ODEs in the Symbolic Math Toolbox documentation for more information about using dsolve.