Ode45 solves an equation that containing a definite integral term

Hi, i have a problem on the following equation when solved by Ode45, which contains a definite integral term. I dont know how to transform it so that it can be solved by ode45.

can someone help me ?

Thank you in advance

[Star Strider]

Is ‘y’ a function of x or t?

[Xuan Ling Zhang]

Thank you for the concern. It is a superscript (') just looks like (t), due to the display problem.

Answer 1

[Torsten]

function main
  a = ...;
  b = ...;
  c = ...; 
  d = ...;
  u0 = 1;
  usol = fzero(@(u)fun(u,a,b,c,d),u0);
  fun_ode = @(t,y)[y(2);y(3);y(4);-usol*y(3)-y(1)^2];
  y0 = [a;b;c;d];
  tspan = [0 1];
  [T,Y] = ode45(fun_ode,tspan,y0);
  plot(T,Y)
end


function res = fun(u,a,b,c,d)
  fun_ode = @(t,y)[y(2);y(3);y(4);-u*y(3)-y(1)^2;y(2)^2];
  y0 = [a;b;c;d;0];
  tspan = [0 1];
  [T,Y] = ode45(fun_ode,tspan,y0);
  res = Y(end,5)-u;
end

[Xuan Ling Zhang]

Dear Torsten，

Thank you for your valuable answer！！！

[sitaram sahu]

Dear Xuan Ling Zhang, please tell how you solved it,

and Dear Torsten, please explain your answer. it would be very helpfull

[Arpan Laha]

Dear Torsten,

Please correct me if I am wrong.

fun_ode = @(t,y)[y(2);y(3);y(4);-u*y(3)-y(1)^2;y(2)^2]; solves the function provided by Xuan replacing the definite integral term by indefinite integral. The y we got as solution from res will not be the same if solved by taking definite integral. Lets term this as Yindf. lets term the actual solution as Ydef. For obvious reason Yindf~=Ydef. Then how can we substitute the value of u into the actual equation ? Thanks

[Torsten]

For a given value of u, you determine the solution y of the differential equation y''''+u*y''+y^2=0 in the interval [0;1] (components 1-4 of fun_ode).

Simultaneously, you integrate y'^2 in the interval [0;1] (component 5 of fun_ode).

Usually, u will not be equal to component 5 of fun_ode, evaluated at x=1.

So, fzero must be used to adjust u such that the two numbers become equal.

Answer 2

[yuan bin]

refer idsolver

https://www.sciencedirect.com/science/article/abs/pii/S0010465513003093?via%3Dihub

IDSOLVER: A general purpose solver for nth-order integro-differential equations

https://cpc.cs.qub.ac.uk/summaries/AEQU_v1_0.html

My code is below:

[x ,nominales] = ode45(@model0,[0,1],[a,b,c,d]);

nominales = [x, nominales];

Tol = 1e-8;

st.TolQuad = 1e-8;

% Iterative solution

error = 1e3;

iteration = 1;

fprintf(' Error convergence\n ');

fprintf(' ================= \n');

fprintf(' Iteration Error \n');

while error > Tol

st.nominales = nominales;

S = ode45(@(x,y)model(x,y,st),[0,1],[a,b,c,d]);

x = nominales(:,1);%time points

R = deval(S, x);

y = R.';

error = sum((y(:,1)-nominales(:,2)).^2);

fprintf(' %4i %8.2e\n',[iteration error]);

alpha = 0.5;

nominales(:,2) = (1-alpha)*nominales(:,2)+alpha*y(:,1);

nominales(:,1) = x;

iteration = iteration+1;

end

warning('on')

plot(x,y);

% legend('y','dy','d2y','d3y');

disp('done');

function dy = model0(x,y)

% Initial guess generator

dy = zeros(4,1);

% y-> [y,dy d2y d3y];

dy(1) = y(2);

dy(2) = y(3);

dy(3) = y(4);

dy(4) = -(y(1)^2+ y(3)* 0); % 0 for initiation value

end

function dy = model(x,y,st)

nominales = st.nominales;

TolQuad = st.TolQuad;

% Interpolation step

ys = @(s) interp1(nominales(:,1),nominales(:,3),s);

% Integro-differential equation

dy = zeros(4,1);

dy(1) = y(2);

dy(2) = y(3);

dy(3) = y(4);

dy(4) = -(y(1)^2 + y(3) *quadl(@(s) ys(s).*ys(s) ,0,1,TolQuad));

end