Ode45 solves an equation that containing a definite integral term

17 次查看(过去 30 天)
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.
Eq.gif
can someone help me ?
Thank you in advance

回答(2 个)

Torsten
Torsten 2019-9-6
编辑:Torsten 2019-9-6
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 个评论
Arpan Laha
Arpan Laha 2021-5-20
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
Torsten 2021-5-20
编辑:Torsten 2021-5-20
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
yuan bin 2023-1-11
编辑:yuan bin 2023-1-11
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

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by