I am trying to solve a coupled set of ODE's different ways. (1) using ode45 which I have been able to do (2) solve them using Euler's method and (3) Heun's method

Question

Dominic Riepenhoff 2022-4-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1699455-i-am-trying-to-solve-a-coupled-set-of-ode-s-different-ways-1-using-ode45-which-i-have-been-able-t

编辑： Torsten 2022-4-19

I am unable to get my Huen's method to properly fit the same like as the ode45 solved for. I do not know how to properly call the variable that was just solved for inside the expression to then resolve for a new value. If this does not make sense i can try to explain further. It follows this rule.

y(i+1)=yi+2(k1+k2)(h)

k_{1}=f(x_{i},y_{i})k1=f(xi,yi)

k_{2}=f(x_{i}+h, y_{i}+k_{1}h)k2=f(xi+h,yi+k1h)

Where for my program y is N (number of moles) and instead of only having one equation i have three coupled equations to solve simaltaneously.

clear
clc
global k1 k2 k3 V0 v0 CA0
k1 = 0.5;
k2 = 0.06;
k3 = 0.04;
V0 = 50;
v0 = 1.5;
CA0 = 0.8;
% ode45 soltion method
tspan = [0:5:200];
N0 = [ 0  0  0 ]; % initial conditions at t=0
[t, N_O] = ode45(@proj2ode45, tspan, N0);
for i = 1:length(t)
    NA(i,1) = N_O(i,1);
    NB(i,1) = N_O(i,2);
    NC(i,1) = N_O(i,3);
    
    
    
end
% Euler Method of solving N
h = 1;
t_E = [0:h:200];
Na = zeros(size(t_E));
Nb = zeros(size(t_E));
Nc = zeros(size(t_E));
N_E =[Na' Nb' Nc'];
n = numel(Na);
for i=1:n-1
    q(i,:) = proj2ode45(t_E(i),N_E(i,:));
    Na(i+1)=Na(i)+h*q(i,1);
    Nb(i+1)=Nb(i)+h*q(i,2);
    Nc(i+1)=Nc(i)+h*q(i,3);
    N_E(i+1,:)=[Na(i+1),Nb(i+1),Nc(i+1)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%source of error%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Heun's Method of solving N
h = .1;
t_H = [0:h:200];
Na_H = zeros(size(t_H));
Nb_H = zeros(size(t_H));
Nc_H = zeros(size(t_H));
Na_s = zeros(size(t_H));
Nb_s = zeros(size(t_H));
Nc_s = zeros(size(t_H));
N_H =[Na_H' Nb_H' Nc_H'];
N_s =[Na_s' Nb_s' Nc_s'];
n = numel(Na_H);
for i=1:n-1
    w(i,:) = proj2ode45(t_H(i),N_H(i,:));
    p(i,:) = proj2ode45(t_H(i),N_H(i,:)); 
    
    % i have tried this method by calling another variable
    Na_s(i+1)=Na_H(i)+h*(w(i,1));
    Nb_s(i+1)=Nb_H(i)+h*(w(i,2));
    Nc_s(i+1)=Nc_H(i)+h*(w(i,3));
    N_s(i+1,:)=[Na_s(i+1),Nb_s(i+1),Nc_s(i+1)];
    
    Na_H(i+1)=Na_H(i)+h*((w(i,1)+(Na_s(i+1)))/2);
    Nb_H(i+1)=Nb_H(i)+h*((w(i,2)+(Nb_s(i+2)))/2);
    Nc_H(i+1)=Nc_H(i)+h*((w(i,3)+(Nc_s(i+3)))/2);
    
    % and i have tried this method by trying to just write in the
    % expression 
    
    Na_H(i+1)=Na_H(i)+h*((w(i,1)+(Na_H(i)+h*(p(i,1))))/2);
    Nb_H(i+1)=Nb_H(i)+h*((w(i,2)+(Nb_H(i)+h*(p(i,2))))/2);
    Nc_H(i+1)=Nc_H(i)+h*((w(i,3)+(Nc_H(i)+h*(p(i,3))))/2);
    
    N_H(i+1,:)=[Na_H(i+1),Nb_H(i+1),Nc_H(i+1)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%^^source of error^^%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ode45
V = V0 + v0*t;
CA = NA./V;
CB = NB./V;
CC = NC./V;
% Euler
V_E = V0 + v0*t_E;
Ca = Na'./V_E';
Cb = Nb'./V_E';
Cc = Nc'./V_E';
% Heun's
V_H = V0 + v0*t_H;
Ca_H = Na_H'./V_H';
Cb_H = Nb_H'./V_H';
Cc_H = Nc_H'./V_H';
[max_CB, j] = max(CB);
t_max = t(j);
help = table(Ca,Cb,Cc);
disp(help);
% ode45 max solution output 
fprintf('The maximum concentration of B is %1.4f mol/L at a time of %2.0f seconds \n\n', max_CB, t_max);
Conc_a = Ca;
Conc_b = Cb;
Conc_c = Cc;
euler = table(t_E',Conc_a,Conc_b,Conc_c);
disp(euler)
conc = table(t,CA, CB, CC);
conc.Properties.VariableNames{'t'} = 'Time (sec)';
conc.Properties.VariableNames{'CA'} = 'Conc of A (mol/L)';
conc.Properties.VariableNames{'CB'} = 'Conc of B (mol/L)';
conc.Properties.VariableNames{'CC'} = 'Conc of C (mol/L)';
disp(conc);
figure(1)
x1 = t;
x2 = t_E;
x3 = t_H;
y1 = CA;
y2 = Ca;
y3 = Ca_H;
plot(x1,y1,'O',x2,y2,'X',x3,y3,'*')
function dNdt = proj2ode45(t,N)
global k1 k2 k3 V0 v0 CA0
NA = N(1);
NB = N(2);
NC = N(3);
V = V0 +v0*t;
CA = NA/V;
CB = NB/V;
dNdt(1) = (CA0*v0)-(V*k1*CA^2)+(V*k2*CB);
dNdt(2) = (V*k1*CA^2)-(V*k2*CB)-(V*k3*CB^2);
dNdt(3) = (V*k3*CB^2);
dNdt = dNdt';
end

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Dominic Riepenhoff 2022-4-19

@Torsten

If you could help again that would be great! You were amazing last time!

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Davide Masiello 2022-4-19

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1699455-i-am-trying-to-solve-a-coupled-set-of-ode-s-different-ways-1-using-ode45-which-i-have-been-able-t#answer_945645

编辑：Davide Masiello 2022-4-19

在 MATLAB Online 中打开

Try the code below

clear,clc

%%%%%%%%%%%%%%%%%%%%%%%%% constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%

k1 = 0.5;

k2 = 0.06;

k3 = 0.04;

V0 = 50;

v0 = 1.5;

CA0 = 0.8;

%%%%%%%%%%%%%%%%%%%%%%%%%% ode45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tspan = [0:5:200];

N0 = [0 0 0];

[to,No] = ode45(@(t,N)proj2ode45(t,N,k1,k2,k3,V0,v0,CA0),tspan,N0);

%%%%%%%%%%%%%%%%%%%%%%%%%% Euler %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

he = 1;

te = [0:he:200];

Ne = zeros(3,length(te));

for i = 1:length(te)-1

q = proj2ode45(te(i),Ne(:,i),k1,k2,k3,V0,v0,CA0);

Ne(:,i+1) = Ne(:,i)+he*q;

end

%%%%%%%%%%%%%%%%%%%%%%%%%% Heunn %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hh = .1;

th = [0:hh:200];

Nh = zeros(3,length(th));

for i = 1:length(th)-1

q1 = proj2ode45(th(i),Nh(:,i),k1,k2,k3,V0,v0,CA0);

Nhbar = Nh(:,i)+hh*q1;

q2 = proj2ode45(th(i+1),Nhbar,k1,k2,k3,V0,v0,CA0);

Nh(:,i+1) = Nh(:,i)+hh/2*(q1+q2);

end

%%%%%%%%%%%%%%%%%%%%%%%%%% Plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)

subplot(1,3,1)

plot(te,Ne(1,:),'oy',to,No(:,1),'k')

title('N_A')

legend('Euler','ode45')

subplot(1,3,2)

plot(te,Ne(2,:),'or',to,No(:,2),'k')

title('N_B')

legend('Euler','ode45')

subplot(1,3,3)

plot(te,Ne(3,:),'og',to,No(:,3),'k')

title('N_C')

legend('Euler','ode45')

figure(2)

subplot(1,3,1)

plot(th,Nh(1,:),'oy',to,No(:,1),'k')

title('N_A')

legend('Heunn','ode45')

subplot(1,3,2)

plot(th,Nh(2,:),'or',to,No(:,2),'k')

title('N_B')

legend('Heunn','ode45')

subplot(1,3,3)

plot(th,Nh(3,:),'og',to,No(:,3),'k')

title('N_C')

legend('Heunn','ode45')

%%%%%%%%%%%%%%%%%%%%%%%%% Function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dNdt = proj2ode45(t,N,k1,k2,k3,V0,v0,CA0)

NA = N(1);

NB = N(2);

NC = N(3);

V = V0 +v0*t;

CA = NA/V;

CB = NB/V;

dNdt(1) = (CA0*v0)-(V*k1*CA^2)+(V*k2*CB);

dNdt(2) = (V*k1*CA^2)-(V*k2*CB)-(V*k3*CB^2);

dNdt(3) = (V*k3*CB^2);

dNdt = dNdt';

end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Torsten 2022-4-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1699455-i-am-trying-to-solve-a-coupled-set-of-ode-s-different-ways-1-using-ode45-which-i-have-been-able-t#answer_945705

编辑：Torsten 2022-4-19

在 MATLAB Online 中打开

global k1 k2 k3 V0 v0 CA0
k1 = 0.5;
k2 = 0.06;
k3 = 0.04;
V0 = 50;
v0 = 1.5;
CA0 = 0.8;
% ode45 soltion method
tspan = [0:5:200];
N0 = [ 0  0  0 ]; % initial conditions at t=0
[t, N_O] = ode45(@proj2ode45, tspan, N0);
for i = 1:length(t)
    NA(i,1) = N_O(i,1);
    NB(i,1) = N_O(i,2);
    NC(i,1) = N_O(i,3);   
end
% Euler Method of solving N
h = 1;
t_E = [0:h:200];
Na = zeros(size(t_E));
Nb = zeros(size(t_E));
Nc = zeros(size(t_E));
N_E =[Na' Nb' Nc'];
n = numel(Na);
for i=1:n-1
    q(i,:) = proj2ode45(t_E(i),N_E(i,:));
    Na(i+1)=Na(i)+h*q(i,1);
    Nb(i+1)=Nb(i)+h*q(i,2);
    Nc(i+1)=Nc(i)+h*q(i,3);
    N_E(i+1,:)=[Na(i+1),Nb(i+1),Nc(i+1)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%source of error%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Heun's Method of solving N
h = .1;
t_H = [0:h:200];
Na_H = zeros(size(t_H));
Nb_H = zeros(size(t_H));
Nc_H = zeros(size(t_H));
%Na_s = zeros(size(t_H));
%Nb_s = zeros(size(t_H));
%Nc_s = zeros(size(t_H));
N_H =[Na_H' Nb_H' Nc_H'];
%N_s =[Na_s' Nb_s' Nc_s'];
n = numel(Na_H);
for i=1:n-1
    w(i,:) = proj2ode45(t_H(i),N_H(i,:));
    %p(i,:) = proj2ode45(t_H(i),N_H(i,:)); 
    
    % i have tried this method by calling another variable
    Na_s = Na_H(i)+h*w(i,1);
    Nb_s = Nb_H(i)+h*w(i,2);
    Nc_s = Nc_H(i)+h*w(i,3);
    N_s = [Na_s,Nb_s,Nc_s];
    %N_s(i+1,:)=[Na_s(i+1),Nb_s(i+1),Nc_s(i+1)];
    
    p(i,:) = proj2ode45(t_H(i+1),N_s);
    Na_H(i+1) = Na_H(i)+h*(w(i,1)+p(i,1))/2;
    Nb_H(i+1) = Nb_H(i)+h*(w(i,2)+p(i,2))/2;
    Nc_H(i+1) = Nc_H(i)+h*(w(i,3)+p(i,3))/2;
    
    % and i have tried this method by trying to just write in the
    % expression 
    
    %Na_H(i+1)=Na_H(i)+h*((w(i,1)+(Na_H(i)+h*(p(i,1))))/2);
    %Nb_H(i+1)=Nb_H(i)+h*((w(i,2)+(Nb_H(i)+h*(p(i,2))))/2);
    %Nc_H(i+1)=Nc_H(i)+h*((w(i,3)+(Nc_H(i)+h*(p(i,3))))/2);
    
    N_H(i+1,:)=[Na_H(i+1),Nb_H(i+1),Nc_H(i+1)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%^^source of error^^%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ode45
V = V0 + v0*t;
CA = NA./V;
CB = NB./V;
CC = NC./V;
% Euler
V_E = V0 + v0*t_E;
Ca = Na'./V_E';
Cb = Nb'./V_E';
Cc = Nc'./V_E';
% Heun's
V_H = V0 + v0*t_H;
Ca_H = Na_H'./V_H';
Cb_H = Nb_H'./V_H';
Cc_H = Nc_H'./V_H';
[max_CB, j] = max(CB);
t_max = t(j);
help = table(Ca,Cb,Cc);
disp(help);
% ode45 max solution output 
fprintf('The maximum concentration of B is %1.4f mol/L at a time of %2.0f seconds \n\n', max_CB, t_max);
Conc_a = Ca;
Conc_b = Cb;
Conc_c = Cc;
euler = table(t_E',Conc_a,Conc_b,Conc_c);
disp(euler)
conc = table(t,CA, CB, CC);
conc.Properties.VariableNames{'t'} = 'Time (sec)';
conc.Properties.VariableNames{'CA'} = 'Conc of A (mol/L)';
conc.Properties.VariableNames{'CB'} = 'Conc of B (mol/L)';
conc.Properties.VariableNames{'CC'} = 'Conc of C (mol/L)';
disp(conc);
figure(1)
x1 = t;
x2 = t_E;
x3 = t_H;
y1 = CA;
y2 = Ca;
y3 = Ca_H;
plot(x1,y1,'O',x2,y2,'X',x3,y3,'*')
function dNdt = proj2ode45(t,N)
global k1 k2 k3 V0 v0 CA0
NA = N(1);
NB = N(2);
NC = N(3);
V = V0 +v0*t;
CA = NA/V;
CB = NB/V;
dNdt(1) = (CA0*v0)-(V*k1*CA^2)+(V*k2*CB);
dNdt(2) = (V*k1*CA^2)-(V*k2*CB)-(V*k3*CB^2);
dNdt(3) = (V*k3*CB^2);
dNdt = dNdt';
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

I am trying to solve a coupled set of ODE's different ways. (1) using ode45 which I have been able to do (2) solve them using Euler's method and (3) Heun's method

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

更多回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

I am trying to solve a coupled set of ODE's different ways. (1) using ode45 which I have been able to do (2) solve them using Euler's method and (3) Heun's method

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

更多回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论