Range-Kutta solving error.

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');

Answer 1

[David Wilson]

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.

[Arkajyoti Chaterjee]

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]

Try increasing the time span such that tfinal is say 1000 and see if you can still integrate it.