Is it possible to solve this differential equation by modifying a given algorithm?
Ältere Kommentare anzeigen
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.
Akzeptierte Antwort
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 Kommentare
Karolis
am 11 Okt. 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
am 11 Okt. 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.
Weitere Antworten (0)
Kategorien
Mehr zu Assumptions finden Sie in Hilfe-Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
