Help in looping and iterating constant values in differential equation which is being solved my numerical method.
1 view (last 30 days)
Show older comments
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')
0 Comments
Answers (1)
Torsten
on 11 May 2022
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 %%%%
15 Comments
Torsten
on 16 May 2022
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));
See Also
Categories
Find more on Geometry and Mesh 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!