System of mixed PDE with 2 variables using FDM and RK4

4 次查看(过去 30 天)
I want solve the attached mixed pdes in 2 variables with FDM and R4K. See my trial below: I received and error that 'N' is uncrecognised.
% parameters.
k1 = 6*10^(-12);
eta1 = 0.0240;
apha3 = -0.001104;
gama1 = 0.1093;
d = 0.0002;
%thetab = 0.0001;
% activity parameter
xi = 0.1;
%%%% Simplify parameters
A = (apha3*k1)/gama1;
B = A/k1;
%%%% Discretize xspace
z_vec = linspace(0, d, N);
dz = z_vec(2) - z_vec(1);
%%%% Discretize t
dt = (dz^2)/(2);
%%%% Discretize t
dtrk4 = (dz^2)/(2);
trk4_vec = 0: dtrk4:10;
% Allocate memory for u
thetark4_mat = zeros(length(z_vec), length(trk4_vec));
vrk4_mat = zeros(length(z_vec), length(trk4_vec));
%%%% IC
thetark4_mat(1,:) = 0; % Left end of pipe
thetark4_mat(end,:) = d; % right end of the pipe
vrk4_mat(1,:) = 0; % Left end of pipe
vrk4_mat(end,:) = d; % right end of the pipe
for tdz = 1:length(trk4_vec)-1
k11 = Derivative(thetark4_mat(:, tdz));
k12 = Derivative(vrk4_mat(:, tdz));
k21 = Derivative(thetark4_mat(:, tdz) + k11*dtrk4/2);
k22 = Derivative(vrk4_mat(:, tdz) + k12*dtrk4/2);
k31 = Derivative(thetark4_mat(:, tdz) + k21*dtrk4/2);
k32 = Derivative(vrk4_mat(:, tdz) + k22*dtrk4/2);
k41 = Derivative(thetark4_mat(:, tdz) + k31*dtrk4);
k42 = Derivative(vrk4_mat(:, tdz) + k32*dtrk4);
phi = (1/6)*(k11+ 2*k21 + 2*k31 +k41);
phi1 = (1/6)*(k12+ 2*k22 + 2*k32 +k42);
thetark4_mat(:, tdz+1) = thetark4_mat(:, tdz) + phi*dtrk4;
vrk4_mat(:, tdz+1) = vrk4_mat(:, tdz) + phi*dtrk4;
end
%%%% Plot this
[tt, zz] = meshgrid(trk4_vec, z_vec);
mesh(zz, tt, thetark4_mat)
xlabel('z')
ylabel('time')
zlabel('theta')
figure()
mesh(zz, tt, vrk4_mat)
xlabel('z')
ylabel('time')
zlabel('v')
function dthetadt = Derivative(theta, v)
%global x_vec dx k
dthetadt = 0*theta;
dvdt = 0*theta;
z_vec = linspace(0,10,10);
% parameters.
k1 = 6*10^(-12);
eta1 = 0.0240;
apha3 = -0.001104;
gama1 = 0.1093;
d = 0.0002;
%thetab = 0.0001;
% activity parameter
xi = 0.1;
%%%% Simplify parameters
A = (apha3*k1)/gama1;
B = A/k1;
dz = z_vec(2) - z_vec(1);
for j = 2:length(z_vec)-1
dvdt(j) = 0;
dthetadt(j) = A/(apha3*dz^2)*(theta(j+1) -2*theta(j) + theta(j-1)) + ...
B/(2*dz)*(v(j+1) - v(j-1));
dvdt(j) = A/(2*dz^3)*(theta(j+2) +2*theta(j+1) -2*theta(j-1)- theta(j-2))+ ...
(eta1 - B)/(dz^2)*(v(j+1) -2*v(j) + v(j-1)) - z/(2*dz)*(theta(j+1) - theta(j-1));
end
end
  20 个评论
Torsten
Torsten 2022-8-6
theta and v both have only 5 elements.
But in these loops
for i=(N+4):(2*N)
rhsode(i,1)= C/(2*h^3)*(v(i+1) + 2*v(i) - 2*v(i-1))+ ...
+ G/(h^3)*(theta(i+2) -2*theta(i+1) + 2*theta(i-1) - theta(i-2))...
+xi/(2*h)*(theta(i+1) - theta(i-1)); %uses central difference for third derivative for theta
end
rhsode(2*N+1,1)=0;
for i = (5+N):(2*N + 1)
rhsode(i,1)= C/(2*h^3)*(v(i+1) + 2*v(i) - 2*v(i-1))+ ...
+ G/(h^3)*(-theta(i-2) +3*theta(i-1) -3*theta(i) + theta(i+1))...
+xi/(2*h)*(theta(i+1) - theta(i-1)); %involves special third derivative
end
almost every index of theta and v exceeds 5.
Use two arrays rhsode_theta and rhsode_v both arrays having 5 elements. Then indexing should be simple.
At the end, return
rhsode = [rhsode_theta;rhsode_v]

请先登录,再进行评论。

回答(2 个)

University Glasgow
% I still don't this error despite having 5 by 5 entries in both you theta and v:
Error in counts_4_IjuptilK_080822 (line 41)
[t,y] = ode15s(@lcode1,tspan,u0, options);
% parameters.
k1 = 6*10^(-12);
eta1 = 0.0240;
apha3 = -0.001104;
gama1 = 0.1093;
d = 0.0002;
N = 4;
% Boundary Conditions
Phi = 0;
% activity parameter
xi = 0.1;
h = d/N; %% step size
%%%% Simplify parameters
A = k1/gama1;
B = apha3/k1;
C = eta1 - B;
G = apha3*A;
Theta = 0.0001;
z=linspace(0,d,N+1);
theta0 = Theta*sin(pi*z/d)
v0 = zeros(1,N+1)
M1=eye(N+1,N+1);
M1(1,1)=0;
M1(N+1,N+1)=0;
M2=zeros(N+1,N+1);
M=[M1 M2;M2 M2];
u0 = [theta0'; v0'];
tspan = [0 10];
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
%[t,y] = ode15s(@lcode1,tspan,[u0;v0], options);
[t,y] = ode15s(@lcode1,tspan,u0, options);
%theta = y(:,1:5);
%v = y(:,1:5);
% Plot the solution.
function rhsode = lcode1(~,y, k1, eta1, apha3, d, Phi, xi, A, B, C, G)
theta = y(1:5);
v = y(6:10);
rhsode = [theta(1)-Phi
(A/(h^2))*(theta(3)+2*theta(2)-theta(1))+(B/(2*h))*(v(3)-v(1))
(A/(h^2))*(theta(4)+2*theta(3)-theta(2))+(B/(2*h))*(v(4)-v(2))
(A/(h^2))*(theta(5)+2*theta(4)-theta(3))+(B/(2*h))*(v(5)-v(3))
theta(5)-Phi
v(N+1)
(G/(h^3))*(-theta(5) + 3*theta(2)-3*theta(3)+theta(4))+(C/(2*h^2))*(v(3) +2*v(2) - v(1)) + (xi/(2*h))*(theta(3)-theta(1))
(G/(h^3))*(-theta(2) + 3*theta(3)-3*theta(4)+theta(5))+(C/(2*h^2))*(v(4) -2*v(3) - v(2)) + (xi/(2*h))*(theta(4)-theta(2))
(G/(h^3))*(-theta(3) + 3*theta(4)-3*theta(5)+theta(2))+(C/(2*h^2))*(v(5) -2*v(4) - v(3)) + (xi/(2*h))*(theta(5)-theta(3))
v(5)];
end
  1 个评论
Torsten
Torsten 2022-8-8
I don't understand why you deleted the call to "fsolve" in an earlier code.
Here is the error again I already commented on.
% parameters.
k1 = 6*10^(-12);
eta1 = 0.0240;
apha3 = -0.001104;
gama1 = 0.1093;
d = 0.0002;
N = 4;
% Boundary Conditions
Phi = 0;
% activity parameter
xi = 0.1;
h = d/N; %% step size
%%%% Simplify parameters
A = k1/gama1;
B = apha3/k1;
C = eta1 - B;
G = apha3*A;
Theta = 0.0001;
z=linspace(0,d,N+1);
theta0 = Theta*sin(pi*z/d)
v0 = zeros(1,N+1)
M1=eye(N+1,N+1);
M1(1,1)=0;
M1(N+1,N+1)=0;
M2=zeros(N+1,N+1);
M=[M1 M2;M2 M2];
u0 = [theta0'; v0'];
tspan = [0 10];
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
%[t,y] = ode15s(@lcode1,tspan,[u0;v0], options);
[t,y] = ode15s(@(t,y)lcode1(t,y,k1, eta1, apha3, d, Phi, h, xi, A, B, C, G),tspan,u0, options);
%theta = y(:,1:5);
%v = y(:,1:5);
% Plot the solution.
function rhsode = lcode1(~,y, k1, eta1, apha3, d, Phi, h, xi, A, B, C, G)
theta = y(1:5);
v = y(6:10);
rhsode = [theta(1)-Phi
(A/(h^2))*(theta(3)+2*theta(2)-theta(1))+(B/(2*h))*(v(3)-v(1))
(A/(h^2))*(theta(4)+2*theta(3)-theta(2))+(B/(2*h))*(v(4)-v(2))
(A/(h^2))*(theta(5)+2*theta(4)-theta(3))+(B/(2*h))*(v(5)-v(3))
theta(5)-Phi
v(1)
(G/(h^3))*(-theta(5) + 3*theta(2)-3*theta(3)+theta(4))+(C/(2*h^2))*(v(3) +2*v(2) - v(1)) + (xi/(2*h))*(theta(3)-theta(1))
(G/(h^3))*(-theta(2) + 3*theta(3)-3*theta(4)+theta(5))+(C/(2*h^2))*(v(4) -2*v(3) - v(2)) + (xi/(2*h))*(theta(4)-theta(2))
(G/(h^3))*(-theta(3) + 3*theta(4)-3*theta(5)+theta(2))+(C/(2*h^2))*(v(5) -2*v(4) - v(3)) + (xi/(2*h))*(theta(5)-theta(3))
v(5)];
end

请先登录,再进行评论。


University Glasgow
Hi, thank you for your help. The code is working well for 4 counts.
But when I tried to generalised, I received this error: Error using daeic12 This DAE appears to be of index greater than 1.
Below is my code:
clear
clc
close all
% parameters.
global k1 eta1 alpha3 gamma1 d N Phi xi h A B C G
k1 = 6*10^(-12);
eta1 = 0.0240;
xi = -0.1; % activity parameter
alpha3 = -0.001104;
gamma1 = 0.1093;
d = 0.0002;
N = 4;
% Boundary Conditions
Phi = 0;
%%%% Simplify parameters
A = k1/gamma1;
B = alpha3/k1;
C = eta1 - B;
G = alpha3*A;
% step size
h = d/N;
% range of z
z=linspace(0,d,N+1);
% initials
Theta = 0.0001;
theta0 = Theta*sin(pi*z/d);
v0 = zeros(1,N+1);
theta0_int=theta0(1:N);
v0_int=v0(1:N);
u0 = [theta0_int'; v0_int'];
% Matrix M
M1=eye(N,N);
M2=zeros(N,N);
M=[M1 M2;M2 M2];
% t span
tspan = [0 4];
% ode solver
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,y] = ode15s(@lcode1,tspan,u0, options);
% Extract the solution for theta and v
%theta = [Phi*ones(length(t), 1) y(:,1:N-1) Phi*ones(length(t), 1)];
%v = [zeros(length(t), 1) y(:,N:2*(N-1)) zeros(length(t), 1)];
% Right hand side of the ODEs: F(U)
function rhsode = lcode1(t,y)
global Phi xi h A B C G N
% initialize theta and v
theta = y(1:N);
v = y(N+1:2*N);
% theta equation
% discretise equation:(A/(h^2))*(theta(i+1)+2*theta(i)-theta(i-1))+(B/(2*h))*(v(i+1)-v(i-1))
for i=2:N
theta(N+1) = 0;
v(N+1) = 0;
rhsode(i,1)=(A/(h^2))*(theta(i+1)+2*theta(i)-theta(i-1))+(B/(2*h))*(v(i+1)-v(i-1)); % internal nodes for 0<z<d
end
% v equation
% for i = N: corresponds to i=9-5 =4 =N
rhsode(N+2,1) = C/(h^2)*(v(N-1) + 2*v(N-2) - v(N-3))+ ...
+ G/(h^3)*(-theta(N-3) +3*theta(N-2) - 3*theta(N-1) + theta(N))...
+xi/(2*h)*(theta(N-1) - theta(N-3)); %involves special third derivative
% for i = 8, this corresponds i = 8-5 =3
% let j =N+4-(N+1) = 3
for j=3:N-1
theta(N+1) = 0;
v(N+1) = 0;
rhsode(j,1)= C/(h^2)*(v(j+1) + 2*v(j) - v(j-1))+ ...
+ G/(2*h^3)*(v(j+2) -2*theta(j+1) + 2*theta(j-1) - theta(j-2))...
+ xi/(2*h)*(theta(j+1) - theta(j-1)); %uses central difference for third derivative for theta
end
% for i = N: corresponds to i=9-5 =4 =N
theta(N+1) = 0;
v(N+1) = 0;
rhsode(2*N,1) = C/(h^2)*(v(N+1) - 2*v(N) + 2*v(N-1))+ ...
+ G/(h^3)*(-theta(N-2) +3*theta(N-1) -3*theta(N) + theta(N+1))...
+xi/(2*h)*(theta(N+1) - theta(N-1)); %involves special third derivative
end

类别

Help CenterFile Exchange 中查找有关 Manage Products 的更多信息

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by