Range-Kutta solving error.

1 次查看（过去 30 天）
Arkajyoti Chaterjee
评论： David Wilson 2021-1-7
I have been asked to solve the set of ODE over (https://imgur.com/a/gjK8H9j) with RK4 using any step-value. Then, compare it with (1) result from any package and (2) the analytical function. But, the RK4/ode45 output seems to be completely botched up.
% Range Kutta Order 4 Code
clc
syms t
y = zeros(1,length(x)); % Pre-allocate array.
y(1) = 10; % Initial Condition.
S = solve((10 - (10+t)*exp(-t)) + 10*exp(-200*t) == y(1), t); % Analytical function solved at initial function so that the start of domain can be found.
h = 0.01; % Step Size - VARY at discretion.
x = S:h:S+1; % Take domain to be 1 (VARY at discretion) and divide domain discretely.
Y = @(r,t) -200*(r - (10 - (10+t)*exp(-t))) + exp(-t)*(9 + t); % Function - VARY at discretion. F has been substituted.
for i=1:(length(x)-1) % Iteration loop. See [https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#The_Runge%E2%80%93Kutta_method] for particular expressions.
k_1 = Y(x(i),y(i));
k_2 = Y(x(i)+0.5*h,y(i)+0.5*h*k_1);
k_3 = Y((x(i)+0.5*h),(y(i)+0.5*h*k_2));
k_4 = Y((x(i)+h),(y(i)+k_3*h));
y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
end
% ODE45
tspan = [0,1]; y0 = 10;
[tx, yx] = ode45(Y, tspan, y0);
% Validation via graph of RK4, ODE45 and Analytical function
plot(x,y,'o-', tx, yx, '--');
fplot(@(t) (10 - (10+t).*exp(-t)) + 10.*exp(-200.*t), [0,1], '-.*c');

请先登录，再回答此问题。

回答（1 个）

David Wilson
David Wilson 2021-1-7
编辑：David Wilson 2021-1-7
Seems fine to me:
Let's try to check the analhytical solution, and also find the derivative of F(t). That's pretty easy:
clc
syms t real
syms y(t)
F = 10-(10+t)*exp(-t);
Fdot = diff(F,t)
ydot = -200*(y-F)+Fdot;
ysoln = dsolve(diff(y,t) == ydot, y(0) == 10) % analytical solution
yana = matlabFunction(ysoln)
Now we can try an analytical solution. We need to know a time span, which you've neglected to give us, so I'm taking t=0 to 10.
I'm using ode45, which is *not* RK45, but close
%% Now try numerical solution
yd = @(t,y) 2000 - 200*y - 199*exp(-t)*(t + 10) - exp(-t)
tspan = [0,10]';
y0 = 10;
[t,y] = ode45(yd,tspan,y0)
plot(t,[y, yana(t)])
Now the real problem is that you have a stiff system, so RK is not ever going to work. (I suspect that is the point of hte homework?)
You will need a stiff system solver such as odexxs. If you get an error using an explicit single-step solver, then that's to be expected.
You can see it is stiff by the parameters in the exponential terms.
  2 个评论
Arkajyoti Chaterjee
Arkajyoti Chaterjee 2021-1-7
Thanks. But, setting h to even lower values tackles stiffness well enough. (The botched output was weirdly corrected after replacing @(r,t) with @(t,r).)
David Wilson
David Wilson 2021-1-7
Try increasing the time span such that tfinal is say 1000 and see if you can still integrate it.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by