# Range-Kutta solving error.

5 views (last 30 days)
Arkajyoti Chaterjee on 7 Jan 2021
Commented: David Wilson on 7 Jan 2021
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');

David Wilson on 7 Jan 2021
Edited: David Wilson on 7 Jan 2021
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 CommentsShowHide 1 older comment
David Wilson on 7 Jan 2021
Try increasing the time span such that tfinal is say 1000 and see if you can still integrate it.

### Categories

Find more on Calculus in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!