Ode45 solves an equation that containing a definite integral term
20 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (2)
Torsten
le 6 Sep 2019
Modifié(e) : Torsten
le 6 Sep 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 commentaires
Arpan Laha
le 20 Mai 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
le 20 Mai 2021
Modifié(e) : Torsten
le 20 Mai 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
le 11 Jan 2023
Modifié(e) : yuan bin
le 11 Jan 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 commentaires
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!