How do I add additional values for k into the function?

clear; clc; close('all');
y = 0.2;
t = 0;
h = 0.01; 
tn = 0.2;
[t, y] = Euler(t, y, h, tn);

function [t, y] = Euler(t0, y0, h, tn)

    t = (t0:h:tn)';
    y = zeros(size(t));
    y(1) = y0;
    for i = 1:1:length(t)-1
        y(i+1) = y(i) + h*f(y(i),t(i));
    end
    fprintf('The concentration of H2O after 0.2 seconds (kf=30):\n')
    disp(y(21))
end

function dydt = f(y,t)
    k1 = 30;
    c = 0.2;
    b = 0.4;
    a = 0.5;

    dydt = k1*b*a;
end

So for example, the output I have now is:

The concentration of H2O after 0.2 seconds (kf=30): 1.4000

How do I add a k2 and a k3 to the function both with values of 20 and 40 and make my output:

The concentration of H2O after 0.2 seconds (kf=30): 1.4000

The concentration of H2O after 0.2 seconds (kf=20): _______

The concentration of H2O after 0.2 seconds (kf=40): _______

[Walter Roberson]

Is the 30 of kf the same as the 30 of k1 in f() ?

[johnmurdock]

Yes, my apologies

[Walter Roberson]

See http://www.mathworks.com/help/matlab/math/parameterizing-functions.html so you can control the k1 value from outside the f function.

Loop over index of hf parameters. Store the results in a 2D array with one of the dimensions being the index of the hf parameter. Do not emit any output until after all of the calculation is done.

Now you can iterate over hf parameter index for each time step, showing the values as appropriate.

[johnmurdock]

Thank you!

[johnmurdock]

I understand controlling the k1 values, however I can't figure out how to impliment it all in the script file.

Answer 1

[Walter Roberson]

y = 0.2;
t = 0;
h = 0.01; 
tn = 0.2;
kf_vals = [30 20 40];
for Kidx = 1 : length(kf_vals)
    kf = kf_vals(Kidx);
    %we assume that all the t values are the same
    [t, y(:,Kidx)] = Euler(t, y, h, tn, kf);
end

for Kidx = 1 : length(kf_vals)
   kf = kf_vals(Kidx);
   fprintf('The concentration of H2O after %g seconds (kf=%g): %g\n', t(end), kf, y(end, Kidx));
 end

function [t, y] = Euler(t0, y0, h, tn, kf)
    t = (t0:h:tn)';
    y = zeros(size(t));
    y(1) = y0;
    for i = 1:1:length(t)-1
        y(i+1) = y(i) + h*f(y(i), t(i), kf);
    end
end

function dydt = f(y, t, kf)
    c = 0.2;
    b = 0.4;
    a = 0.5;
    dydt = kf*b*a;
end