Is it possible to solve this differential equation by modifying a given algorithm?
Afficher commentaires plus anciens
Hello everyone,
Sorry for bothering with a similar problem (https://www.mathworks.com/matlabcentral/answers/2031639-is-there-a-possibility-to-solve-this-equation-numerically-using-matlab?s_tid=srchtitle) which was previously resolved by a user named Torsten. I am now trying to solve this final version of differential equation using MATLAB (r and u are variables):
The code used for solving this equation:
u0 = 0;
u0dot = 0.001;
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))
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)])
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) 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
It seems to me that the equation is stiff so with each iteration I am constantly getting this error: Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is... The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy. Maybe it is impossible to solve this equation numerically using MATLAB? And once again, thank you for your answers.
Réponse acceptée
Seems your function to be integrated has a singularity between d/2 and R. Applying "integral" in this case is almost impossible. It has nothing to do with "stiffness".
Ecm = 33787;
d = 0.02;
R = 0.085;
Ac = 0.02;
S = 54;
C = 750;
Gcm = 14000;
fun = @(r,y3) y3*S*Ecm./r/(8*(S+C)).*(4*Ac/pi-4*r.^2+d^2)./(Gcm-2223.5*y3*S*Ecm./r/(8*(S+C)).*(4*Ac/pi-4*r.^2+d^2));
y3 = 0:0.1:2;
hold on
for i=1:numel(y3)
plot(linspace(d/2,R,1000),fun(linspace(d/2,R,1000),y3(i)))
end
hold off
4 commentaires
Karolis
le 11 Oct 2023
Yes, the location of the singularity depends on the value of u''(x) (and maybe the singularity in the r-interval disappears if u''(x) is different from 0 <= u''(x) <= 2 ; I didn't check it in detail).
You could theoretically find the zeros of the denominator as a function of u''(x) and see what values for u''(x) are allowed for that the zeros are outside the interval [d/2 , R]. But I don't know if this helps in the solution of the above problem because u'' follows from u and u' - you cannot prescribe it, of course.
Karolis
le 11 Oct 2023
I don't understand what your code from above is good for to solve the original problem.
You have to apply it the opposite way to use it in your problem: Solve f(d2u) - y(1) = 0 for d2u. But I think it doesn't matter much whether you define f as the function above or if you directly use the integral representation.
Plus de réponses (0)
Catégories
En savoir plus sur Assumptions dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
