Solving system of equations

8 visualizaciones (últimos 30 días)
EldaEbrithil
EldaEbrithil el 27 de Mayo de 2020
Comentada: darova el 30 de Mayo de 2020
Hi all
i have a question about solving this system of equations. Tt, Pt and M are related to space and time due to i and j; i want to solve the system maintaining that dependence, so the result will be a matrix respectively for Tt, Pt and M. When i try to solve, i obtain "Out of range subscript." error. gamma, deltax and deltat are constant
Thanks to all
Tt=zeros(length(x),length(t));
Pt=zeros(length(x),length(t));
M=zeros(length(x),length(t));
Tt(1,1)=3.000555630247608e+02;
Pt(1,1)=2.201018491400215e+05;
M(1,1)=0.023565919700319;
for j=1:length(t)-1
for i=2:length(x)-1
Alla = cell(length(x),length(t));
Allb = cell(length(x),length(t));
Allc = cell(length(x),length(t));
syms Tt Pt M
[sola,solb,solc]=vpasolve(Tt(i,j+1)==0.5*(Tt(i+1,j)-Tt(i-1,j))+((1+((gamma-1)/2)*M(i,j)^2)^(gamma/(gamma-1)))*((Tt(i+1,j)-Tt(i-1,j))*deltat/(2*deltax))+((1+((gamma-1)/2)*M(i,j)^2))*((Pt(i+1,j)-Pt(i-1,j))*deltat/(2*deltax)),...
Pt(i,j+1)==0.5*(Pt(i+1,j)-Pt(i-1,j))+2*((1+((gamma-1)/2)*M(i,j)^2)^(gamma/(gamma-1)))*((Tt(i+1,j)-Tt(i-1,j))*deltat/(2*deltax))+3*((1+((gamma-1)/2)*M(i,j)^2))*((Pt(i+1,j)-Pt(i-1,j))*deltat/(2*deltax)),...
M(i,j+1)==0.5*(M(i+1,j)-M(i-1,j))+2*((1+((gamma-1)/2)*M(i,j)^2)^(gamma/(gamma-1)))*((Tt(i+1,j)-Tt(i-1,j))*deltat/(2*deltax))+3*((1+((gamma-1)/2)*M(i,j)^2))*((Pt(i+1,j)-Pt(i-1,j))*deltat/(2*deltax)));
Alla{i,j} = sola;
Allb{i,j} = solb;
Allc{i,j} = solc;
end
end
  17 comentarios
darova
darova el 27 de Mayo de 2020
I can't explain it here
can be re-written as (P(i,j+1)-P(i,j))/dt
can be re-written as (P(i+1,j)-P(i,j))/dx
you what i mean?
Read about this method. Read about "Method of lines"
EldaEbrithil
EldaEbrithil el 28 de Mayo de 2020
Yes i understand, but i think it is what similar to what i have done in my code, the only difference is related to the typology of discretization: you have used a forward discretiation in space and time, i have used a Forward Time Centered Space, FTCS discretization. Thi is the only difference, but the problem i have is easier than you think: i do not understand how to write the code for solving the system of equations practically.

Iniciar sesión para comentar.

Respuestas (1)

darova
darova el 28 de Mayo de 2020
Here is a simple example. I hope it's clear enough. TR, TL, TD - boundary conditions (right, left and down boundaries)
  2 comentarios
EldaEbrithil
EldaEbrithil el 30 de Mayo de 2020
I have tried to implement the method for the equation tht you give me in the.m file but i am not very confident about the results
clc,clear
% problem definition and discretization
dx = 0.01;
dt = 0.008;
xdomain = [0 1];
tdomain = [0 1];
nx = round((xdomain(2)-xdomain(1))/dx);
nt = round((tdomain(2)-tdomain(1))/dt);
x = linspace(xdomain(1),xdomain(2),nx);
t = linspace(tdomain(1),tdomain(2),nt);
u = zeros(nt,nx);
% du/dt - 2*t*du/dx = 0
u(1,:) = sin(2*pi*x);
for k = 1:nt-1
for i = 1:nx-1
% Predictor step
u(k+1,i) = 2*t(k)*dt/dx*(u(k,i+1)-u(k,i)) + u(k,i);
end
end
figure(1);set(gcf,'Visible', 'off')
plot(x,u(85,:))
figure(4);set(gcf,'Visible', 'off')
surf(x,t,u)
%%%%%LAX WENDROFF%%%%%
dx2 = 0.01;
dt2 = 8e-4;
xdomain2 = [0 1];
tdomain2 = [0 1];
nx2 = round((xdomain2(2)-xdomain2(1))/dx2);
nt2 = round((tdomain2(2)-tdomain2(1))/dt2);
x2= linspace(xdomain2(1),xdomain2(2),nx2);
t2 = linspace(tdomain2(1),tdomain2(2),nt2);
u2 = zeros(nx2,nt2);
u2(:,1) = sin(2*pi*x2);%initial condition
for i=2:nx2-1
for j=1:nt2-1
u2(i,j+1)=u2(i,j)+(2*t2(j)*dt2/(2*dx2))*(u2(i+1,j)-u2(i-1,j))+((dt2^2)/(2*dx2))*(u2(i+1,j)-u2(i-1,j))+2*((dt2^2)/(dx2^2))*(t2(j)^2)*(u2(i+1,j)-2*u2(i,j)+u2(i-1,j));
stab(j)=2*t2(j)*dt2/(2*dx2);%always less then one
end
end
figure(2);set(gcf,'Visible', 'off')
plot(x2,u2(:,85))
figure(3);set(gcf,'Visible', 'off')
surf(t2,x2,u2)
darova
darova el 30 de Mayo de 2020
I can't check it. It's too complicated, sorry

Iniciar sesión para comentar.

Categorías

Más información sobre Fluid Dynamics en Help Center y File Exchange.

Productos


Versión

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by