Range-Kutta solving error.

6 views (last 30 days)
Arkajyoti Chaterjee
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');

Answers (1)

David Wilson
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 Comments
Arkajyoti Chaterjee
Arkajyoti Chaterjee on 7 Jan 2021
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 on 7 Jan 2021
Try increasing the time span such that tfinal is say 1000 and see if you can still integrate it.

Sign in to comment.

Categories

Find more on Function Creation in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!