Solving system of equations

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

[darova]

• "Out of range subscript." error. gamma, deltax and deltat are constant

I don't see these variable are defined. Did you define them somewhere?

[EldaEbrithil]

yes

Ltantfinal=0.852059141962213;
Nc=0.1;% ampiezza intervallo 
x=[0:Nc:Ltantfinal];%valori dell'ascissa curvilinea di camera
t=[0:Nc:60];
deltax=x(2)-x(1);
deltat=t(2)-t(1);

[darova]

What is going on here?

• are variable Tt Pt M numerical or symbolic?

[EldaEbrithil]

Symbolic but i need to define Tt Pt and M by zeros due to the for loop

[EldaEbrithil]

Basically i want to solve the system for each i and j

[darova]

I'm confused. What variable are changing and variable you want to find?

[EldaEbrithil]

The variable that i want to find and those that change are both Tt, Pt and M. Is it possible to find these 3 unknowns, considering that are connected in that way thanks to the foor loop? I need to maintain that connection scheme between the unknowns. Are there alternative ways to solve that problem (find Tt, Pt and M, that i repeat, are connected each other) without involving a system solver?

[darova]

Can you show original equations? If you have 3 equations you can get only 3 solutions

[EldaEbrithil]

yes i know that, but i want that to obtain 3 solution for each value of i and j.

[darova]

You have system of PDE (partial differential equations). Do know what does it mean?

[EldaEbrithil]

yes i know, the only way for solving this system is to applied the method of Lax-Wendroff or other reduction methods because i have hyperbolic equations and i can not solve this kind of system throught matlab functions

[darova]

What about FDM(finite difference method)?

[EldaEbrithil]

I have choosen Lax-Wendroff due to his high accuracy; the problem is related especially to the method implementation, i really don't know if the code i wrote is valid for the solution of the system

[darova]

YOu should read about FDM

wikipedia

See simple example attached

[EldaEbrithil]

Ok this is interesting but what about a system, instead of single equation, of equations of this type?

u(k+1,i) = 2*t(k)*dt/dx*(u(k,i+1)-u(k,i)) + u(k,i);

[darova]

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]

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.

Answer 1

[darova]

Here is a simple example. I hope it's clear enough. TR, TL, TD - boundary conditions (right, left and down boundaries)

[EldaEbrithil]

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]

I can't check it. It's too complicated, sorry