Is it possible to solve this differential equation by modifying a given algorithm?
이전 댓글 표시
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.
채택된 답변
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
Karolis
2023년 10월 11일
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
2023년 10월 11일
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.
추가 답변 (0개)
카테고리
도움말 센터 및 File Exchange에서 Assumptions에 대해 자세히 알아보기
Translated by 
