Ode45 solves an equation that containing a definite integral term
18 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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
2 comentarios
Respuestas (2)
Torsten
el 6 de Sept. de 2019
Editada: Torsten
el 6 de Sept. de 2019
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
4 comentarios
Arpan Laha
el 20 de Mayo de 2021
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
el 20 de Mayo de 2021
Editada: Torsten
el 20 de Mayo de 2021
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.
yuan bin
el 11 de En. de 2023
Editada: yuan bin
el 11 de En. de 2023
refer idsolver
IDSOLVER: A general purpose solver for nth-order integro-differential equations
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
0 comentarios
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!