Help in looping and iterating constant values in differential equation which is being solved my numerical method.

Hello guys!! I am solving a differential equation which is time dependent with numerial method but the constants change at every instant of time. In the below code i have kept same constants throughout. Does anyone can help me to write a code to iretrate those constants every time.? I am guessing i need to put lot of for loops.

clear all
% Adams-Bashforth Predictor Corrector Method
% CONSTANTS WHICH VARY WITH TIME
Densityg=5.80;    
Densityl=61.3;
ml=0.1;
T = 21;              
P = 20; 
dorogdoPT=3.40;
dorogdoTP=-2.16;
doroldoPT=18.323;
doroldoTP=-0.51;

%%%% t= mg y =t %%%%

f = @(t,y) (((t/Densityg^2)*((dorogdoPT)*(P/y)+(dorogdoTP)*(T/y)))+((ml/Densityl^2)*((doroldoPT)*(P/y)+(doroldoTP)*(T/y))))/((1/Densityg)-(1/Densityl));     %%%   DIFFERENTIAL EQUATION   %%%%

a = input('Intial value of time (t), a:  ');
b = input('Enter value of time (t) at which mg is to be calculated b:  ');
n = input('Enter no. of subintervals, n: ');
%h = input('Enter the step size,h: ');
alpha = input('Enter the initial condition of mg, alpha:  ');

h = (b-a)/n;
t(1) = a;
w(1) = alpha;
fprintf('    t         w\n');
fprintf('%5.4f  %11.8f\n', t(1), w(1));

for i = 1:3 
    t(i+1) = t(i)+h;
    k1 = h*f(t(i), w(i));
    k2 = h*f(t(i)+0.5*h, w(i)+0.5*k1);
    k3 = h*f(t(i)+0.5*h, w(i)+0.5*k2);
    k4 = h*f(t(i+1), w(i)+k3);
    w(i+1) = w(i)+(k1+2.0*(k2+k3)+k4)/6.0;
    fprintf('%5.4f  %11.8f\n', t(i+1), w(i+1));
end

for i = 4:n 
    t0 = a+i*h;
    part1 = 55.0*f(t(4),w(4))-59.0*f(t(3),w(3))+37.0*f(t(2),w(2));
    part2 = -9.0*f(t(1),w(1));
    w0 = w(4)+h*(part1+part2)/24.0;
    part1 = 9.0*f(t0,w0)+19.0*f(t(4),w(4))-5.0*f(t(3),w(3))+f(t(2),w(2));
    w0 = w(4)+h*(part1)/24.0;
    fprintf('%5.4f  %11.8f\n', t0, w0);
    for j = 1:3 
        t(j) = t(j+1);
        w(j) = w(j+1);
    end
    t(4) = t0;
    w(4) = w0;
end

plot(t,w)
xlabel('time')
ylabel('mg')

Answer 1

[Torsten]

Define all time-dependent parameters as function handles and apply these handles in your function definition.

As an example for Densityg:

Densityg = @(t) 5.8*t;
f = @(t,y) (((t/Densityg(t)^2)*((dorogdoPT)*(P/y)+(dorogdoTP)*(T/y)))+((ml/Densityl^2)*...
           ((doroldoPT)*(P/y)+(doroldoTP)*(T/y))))/((1/Densityg(t))-(1/Densityl));     %%%   DIFFERENTIAL EQUATION   %%%%

[Omkar Gurav]

one more things i would like to add. i have all the values of Densityg with respect to time.

t(s) Densityg

0 5.8

1 5.9

2 7.9

3 6.5 etc

in the same manner i have all the values of all the constants with resp. to time. I need to add these values as the code runs for every second

[Torsten]

Densityg = @(t) interp1(time,densityg,t)

where time, densityg are arrays with

time = tstart:1:tend
densityg(i) = gas density at t=time(i)

[Omkar Gurav]

Hello torsten!

I am trying that but i am not to run the code it getting stuck in the numerical method code line

k1 = h*f(t(i), w(i)); or i guess i am not able to integrate it properly.

Could you please help me in this to itregrate in my code?. Thanks for your time.

[Torsten]

What's your actual code after the changes you made ?

[Omkar Gurav]

clear all

% Adams-Bashforth Predictor Corrector Method

% Approximate the solution to the initial-value problem

% CONSTANTS WHICH VARY WITH TIME

%Densityg=5.80;

densityg = [5.8 5.9 7.9 6.5];

Densityg = @(t) interp1(time,densityg,t)

time = tstart:1:tend

densityg(i) = gasdensity at t=time(i)

%Densityg = 5.8;

Densityl=61.3;

ml=0.1;

T = 21;

P = 20;

dorogdoPT=3.40;

dorogdoTP=-2.16;

doroldoPT=18.323;

doroldoTP=-0.51;

%%%% t= mg y =t %%%%

f = @(t,y) (((t/Densityg^2)*((dorogdoPT)*(P/y)+(dorogdoTP)*(T/y)))+((ml/Densityl^2)*((doroldoPT)*(P/y)+(doroldoTP)*(T/y))))/((1/Densityg)-(1/Densityl));

a = input('Intial value of time (t), a: ');

b = input('Enter value of time (t) at which mg is to be calculated b: ');

n = input('Enter no. of subintervals, n: ');

%h = input('Enter the step size,h: ');

alpha = input('Enter the initial condition of mg, alpha: ');

h = (b-a)/n;

t(1) = a;

w(1) = alpha;

fprintf(' t w\n');

fprintf('%5.4f %11.8f\n', t(1), w(1));

for i = 1:3

t(i+1) = t(i)+h;

k1 = h*f(t(i), w(i));

k2 = h*f(t(i)+0.5*h, w(i)+0.5*k1);

k3 = h*f(t(i)+0.5*h, w(i)+0.5*k2);

k4 = h*f(t(i+1), w(i)+k3);

w(i+1) = w(i)+(k1+2.0*(k2+k3)+k4)/6.0;

fprintf('%5.4f %11.8f\n', t(i+1), w(i+1));

end

for i = 4:n

t0 = a+i*h;

part1 = 55.0*f(t(4),w(4))-59.0*f(t(3),w(3))+37.0*f(t(2),w(2));

part2 = -9.0*f(t(1),w(1));

w0 = w(4)+h*(part1+part2)/24.0;

part1 = 9.0*f(t0,w0)+19.0*f(t(4),w(4))-5.0*f(t(3),w(3))+f(t(2),w(2));

w0 = w(4)+h*(part1)/24.0;

fprintf('%5.4f %11.8f\n', t0, w0);

for j = 1:3

t(j) = t(j+1);

w(j) = w(j+1);

end

t(4) = t0;

w(4) = w0;

end

plot(t,w)

xlabel('time')

ylabel('mg')

[Omkar Gurav]

@Torsten the above is my code. i am beginner in matlab. could you please help in running that code. Thanks again

[Torsten]

The below code modifies the prescription of Densityg.

The result seems questionable since the plot does not start at t=1.

clear all
% Adams-Bashforth Predictor Corrector Method
% Approximate the solution to the initial-value problem
% CONSTANTS WHICH VARY WITH TIME
%Densityg=5.80;
time=1:4;
densityg = [5.8 5.9 7.9 6.5];
Densityg = @(t) interp1(time,densityg,t)



Densityl=61.3;
ml=0.1;
T = 21;
P = 20;
dorogdoPT=3.40;
dorogdoTP=-2.16;
doroldoPT=18.323;
doroldoTP=-0.51;
%%%% t= mg y =t %%%%

f = @(t,y) (((t/Densityg(t)^2)*((dorogdoPT)*(P/y)+(dorogdoTP)*(T/y)))+((ml/Densityl^2)*((doroldoPT)*(P/y)+(doroldoTP)*(T/y))))/((1/Densityg(t))-(1/Densityl));
a = 1;%input('Intial value of time (t), a:  ');
b = 4;%input('Enter value of time (t) at which mg is to be calculated b:  ');
n = 12;%input('Enter no. of subintervals, n: ');
%h = input('Enter the step size,h: ');
alpha = 1.0;%input('Enter the initial condition of mg, alpha:  ');

h = (b-a)/n;
t(1) = a;
w(1) = alpha;
fprintf('    t         w\n');
fprintf('%5.4f  %11.8f\n', t(1), w(1));
for i = 1:3
    t(i+1) = t(i)+h;
    k1 = h*f(t(i), w(i));
    k2 = h*f(t(i)+0.5*h, w(i)+0.5*k1);
    k3 = h*f(t(i)+0.5*h, w(i)+0.5*k2);
    k4 = h*f(t(i+1), w(i)+k3);
    w(i+1) = w(i)+(k1+2.0*(k2+k3)+k4)/6.0;
    fprintf('%5.4f  %11.8f\n', t(i+1), w(i+1));
end
for i = 4:n
    t0 = a+i*h;
    part1 = 55.0*f(t(4),w(4))-59.0*f(t(3),w(3))+37.0*f(t(2),w(2));
    part2 = -9.0*f(t(1),w(1));
    w0 = w(4)+h*(part1+part2)/24.0;
    part1 = 9.0*f(t0,w0)+19.0*f(t(4),w(4))-5.0*f(t(3),w(3))+f(t(2),w(2));
    w0 = w(4)+h*(part1)/24.0;
    fprintf('%5.4f  %11.8f\n', t0, w0);
    for j = 1:3
        t(j) = t(j+1);
        w(j) = w(j+1);
    end
    t(4) = t0;
    w(4) = w0;
end

plot(t,w)
xlabel('time')
ylabel('mg')

[Omkar Gurav]

@Torsten hello! I have a doubt why you used interpl? we dont have any missing data right.

[Torsten]

If you don't choose n = 3, you have missing data.

[Omkar Gurav]

@Torsten hello!. Thanks for your time again.

Need to take n = 4 not 12

n- which is no.of subintervals i am keeping same as b-value of time (t) at which mg is to be calculated b:

i.e. when a =1 b=4 i am keeping n = 4 because i know values at every second keeping n =4 makes step size as 1.

hope the problem statement is more clear now and you could help me.

[Torsten]

If you ask for densityg exactly at the times when the data are supplied (as you claim you will do), you'll get back the exact value for the variable at that time. So there is no reason not to use interp1. Since you only have t and y available when you want to evaluate densityg in your function f, this is the easiest solution.

[Omkar Gurav]

okay @Torsten thanks for that i am trying it with all the constants. But dont know why the results doesnot start at t = 1.

[Torsten]

The Adams part of your code looks wrong to me. This could be the reason for the wrong initial time, especially the part

for j = 1:3 
    t(j) = t(j+1);
    w(j) = w(j+1);
end
t(4) = t0;
w(4) = w0;

[Omkar Gurav]

@Torsten one more thing i am not able to add all time-dependent parameters as function handles and run it. Below i have added one more time dependant variable 'T' but i am getting the following error

t w

0.0000 1.00000000

Error using / Arguments must be numeric, char, or logical.

Error in new_t_code>@(t,y)(((t/Densityg(t)^2)*((dorogdoPT)*(P/y)+(dorogdoTP)*(T/y)))+((ml/Densityl^2)*((doroldoPT)*(P/y)+(doroldoTP)*(T/y))))/((1/Densityg(t))-(1/Densityl))

(line 22)

f = @(t,y)

(((t/Densityg(t)^2)*((dorogdoPT)*(P/y)+(dorogdoTP)*(T/y)))+((ml/Densityl^2)*((doroldoPT)*(P/y)+(doroldoTP)*(T/y))))/((1/Densityg(t))-(1/Densityl));

Error in new_t_code (line 35)

k1 = h*f(t(i), w(i));

clear all
clc
% Adams-Bashforth Predictor Corrector Method
% Approximate the solution to the initial-value problem
% CONSTANTS WHICH VARY WITH TIME
%Densityg=5.80;
time=0:6;
densityg = [5.8 5.9 7.9 6.5 8 9 5];
Densityg = @(t) interp1(time,densityg,t)

Densityl=61.3;
ml=0.1;
T = [11 12 13 14 15 16 17 ];
T = @(t) interpl(time,T,t)
%T = 21;
P = 20;
dorogdoPT=3.40;
dorogdoTP=-2.16;
doroldoPT=18.323;
doroldoTP=-0.51;
%%%% t= mg y =t %%%%
f = @(t,y) (((t/Densityg(t)^2)*((dorogdoPT)*(P/y)+(dorogdoTP)*(T/y)))+((ml/Densityl^2)*((doroldoPT)*(P/y)+(doroldoTP)*(T/y))))/((1/Densityg(t))-(1/Densityl));
a = 0;%input('Intial value of time (t), a:  ');
b = 6;%input('Enter value of time (t) at which mg is to be calculated b:  ');
n = 6;%input('Enter no. of subintervals, n: ');
%h = input('Enter the step size,h: ');
alpha = 1.0;%input('Enter the initial condition of mg, alpha:  ');
h = (b-a)/n;
t(1) = a;
w(1) = alpha;
fprintf('    t         w\n');
fprintf('%5.4f  %11.8f\n', t(1), w(1));
for i = 1:3
    t(i+1) = t(i)+h;
    k1 = h*f(t(i), w(i));
    k2 = h*f(t(i)+0.5*h, w(i)+0.5*k1);
    k3 = h*f(t(i)+0.5*h, w(i)+0.5*k2);
    k4 = h*f(t(i+1), w(i)+k3);
    w(i+1) = w(i)+(k1+2.0*(k2+k3)+k4)/6.0;
    fprintf('%5.4f  %11.8f\n', t(i+1), w(i+1));
end
for i = 4:n
    t0 = a+i*h;
    part1 = 55.0*f(t(4),w(4))-59.0*f(t(3),w(3))+37.0*f(t(2),w(2));
    part2 = -9.0*f(t(1),w(1));
    w0 = w(4)+h*(part1+part2)/24.0;
    part1 = 9.0*f(t0,w0)+19.0*f(t(4),w(4))-5.0*f(t(3),w(3))+f(t(2),w(2));
    w0 = w(4)+h*(part1)/24.0;
    fprintf('%5.4f  %11.8f\n', t0, w0);
    for j = 1:3
        t(j) = t(j+1);
        w(j) = w(j+1);
    end
    t(4) = t0;
    w(4) = w0;
    
plot(t,w); hold on
xlabel('time')
ylabel('mg')
end

[Torsten]

a) Name the function handle different from the name of the array :

Temp = [11 12 13 14 15 16 17 ];
T = @(t) interpl(time,Temp,t);

b) Refer to T as T(t) in your function definition and not as T:

f = @(t,y) (((t/Densityg(t)^2)*((dorogdoPT)*(P/y)+(dorogdoTP)*(T(t)/y)))+((ml/Densityl^2)*((doroldoPT)*(P/y)+(doroldoTP)*(T(t)/y))))/((1/Densityg(t))-(1/Densityl));