Is there a possibility to solve this equation numerically using MATLAB?

Hello,

Is there a possibility to solve a given equation (model) numerically using MATLAB? (u and r are variables).

I have tried solving it by firstly integrating using Symbolic Math Toolbox with respect to r and then using ode15i when x values are from 0 to 0.12. To some extent, this approach works and I get reasonable results. However with a certain initial conditions of u'(0) (typical values of u'(0) are from 0 to 0.005 and u(0) is always equal to 0) the problem becomes stiff and I get this error message: Unable to meet integration tolerances without reducing the step size below the smallest value allowed... Is there something that can be done to solve this problem for a wider variety of initial conditions of u'(0)? When u(x) values (those which were solved by ode15i) were compared to experimental results the model showed a very nice correlation. Maybe I should use a different solver or maybe this problem is unsolvable and I should reconsider my model? Thank you for all advices!

Answer 1

[Torsten]

Please include the code you are using at the moment and critical values for u'(0).

Is there any equation for u without the integral term ? It seems to me that the equation from above is already deduced from some "original" form.

u0 = 0;
u0dot = 0.00001;
u0ddot = 1;  % Estimate
t0 = 0;
y0 = [u0,u0dot,u0ddot];
s = @(x)Ode(t0,[y0(1:2),x]);
h = @(x) subsref(s(x), struct('type', '()', 'subs', {{3}}));
sol = fsolve(h,y0(3),optimset('TolFun',1e-10,'TolX',1e-10))

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

sol = 2.6442e-05

y0 = [y0(1:2),sol];
tspan = [t0,1];
OdeFcn = @Ode;
Mass = [1 0 0;0 1 0;0 0 0];
options = odeset('RelTol',1e-10,'AbsTol',1e-10,'Mass',Mass,'InitialStep',1e-10);
[X,Y] = ode15s(@Ode,tspan,y0,options);
figure(1)
plot(X,[Y(:,1),Y(:,2),Y(:,3)])

for i = 1:numel(X)
  dydx = Ode(X(i),Y(i,:));
  res3(i) = dydx(3);
end
figure(2)
plot(X,res3.')

function dydx = Ode(x,y)
  persistent iter
  if isempty(iter)
    iter = 0;
  end
  iter = iter + 1;
  if mod(iter,10000)==0
    iter = 0;
    x
  end
  Ecm = 33787;
  d = 0.02;
  R = 0.085;
  Ac = 0.02;
  S = 54;
  C = 750;
  Gcm = 14000;

  dydx = zeros(3,1);
  dydx(1) = y(2);
  dydx(2) = y(3);
  fun = @(r)2223.5*y(3)*S*Ecm*r/(8*(S+C)).*(4*Ac/pi-4*r.^2+d^2)./(Gcm-2223.5*y(3)*S*Ecm*r/(8*(S+C)).*(4*Ac/pi-4*r.^2+d^2));
  dydx(3) = y(1) - integral(fun,d/2,R,'AbsTol',1e-10,'RelTol',1e-10);
end

[Karolis]

Wow! Thank you very much Torsten! You saved me tons of time searching for the solution! It works perfectly :) . Now I can fully compare my model to the experimental results.