Solving system of equations

1 次查看（过去 30 天）
EldaEbrithil
EldaEbrithil 2020-5-27
评论： darova 2020-5-30
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 个评论
显示 15更早的评论
darova
darova 2020-5-27
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 2020-5-28
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.

请先登录，再进行评论。

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

回答（1 个）

darova
darova 2020-5-28
Here is a simple example. I hope it's clear enough. TR, TL, TD - boundary conditions (right, left and down boundaries)
  2 个评论
EldaEbrithil
EldaEbrithil 2020-5-30
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 2020-5-30
I can't check it. It's too complicated, sorry

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 General Physics 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by