Help me please with this code:
Afficher commentaires plus anciens
%function [u_save,x,t_save]=sol_direc( M, t_f, tsteps,METODO)
%function [u_save,x,t_save]=sol_direc( M, t_f, tsteps)
M =20;
t_f = 1;
tsteps =1000;
%NMAX = 100;
%resuelve la ecuacion del calor u_t = (au_x)_x en [0,1] con dato inicial
%sin(pi x)
%y condiciones de contorno tipo Dirichlet u(0,t)=u(1,t)=0
%
%Datos de entrada:
% M: numero de partes iguales en que se descompone [0,1]
% t_f: instante de tiempo hasta el que calculamos la solución
% tsteps: numero de partes iguales en que se descompone [0,t_f]
%METODO entre 1 para superficie o 2 para algunos valores
%N=M;
h = 0.5/M;
%Mallado
dt=t_f/tsteps;
tiempos=0:dt:t_f;
%
%h=1/M;
x=0:h:1;
%
% Datos Teóricos:
syms xx tt
us = exp(tt).*sin(pi.*xx);
as = 40.*xx.^2./(40.*xx.^2+1)+1;
fs = diff(us,tt)-diff(as*diff(us,xx),xx);
ai = double(subs(as,xx,x))'; % a numérico
%plot(ai)
% Datos numéricos exactos fi y u
[tm,xm] = meshgrid(tiempos,x);
fi = double(subs(fs,{xx,tt},{xm,tm}));
un=double(subs(us,{xx,tt},{xm,tm}));
%figure(1)
%mesh(un)
%Datos iniciales
xred=x(2:2*M); %quitamos los extremos del vector x
u=sin(pi*xred)'; %trasponemos para obtener un vector columna
%
%
%El operador D
%Usamos el comando diag(vector,k) para crear una matriz tridiagonal
% Matriz tridiagonal,
% b es la diagonal principal, a es la inferior y c es la superior
a =zeros(1,2*M-2);
b =zeros(1,2*M-1);
c =zeros(1,2*M-2);
%
b(1) = (ai(1)+ai(2))/h;
c(1) = -(0.25*ai(1)+0.75*ai(2))/h;
a(1)=-(0.25*ai(1)+0.25*ai(2+1))/h;
a(2*M-2) = -(0.75*ai(M-2)+0.25*ai(M+1))/h;
b(2*M-1) = (ai(M)+ai(M+1))/h;
for i=2:2*M-2
if mod(i,2)== 1
a(i) = -(0.75*ai((i-1)/2)+0.25*ai((i+1)/2))/h;
b(i) = (ai((i-1)/2)+ai((i+1)/2))/h;
c(i) = -(0.25*ai((i-1)/2)+0.75*ai((i+1)/2))/h;
else
a(i) = -(0.25*ai((i-2)/2+1)+0.25*ai(i/2+1))/h;
b(i) = (0.25*ai((i-2)/2+1)+1.5*ai(i/2+1)/h +0.25*ai((i+2)/2+1))/h;
c(i) = -(0.75*ai(i/2+1)+0.25*ai((i+2)/2+1))/h;
end
end
tridiag = diag(a,-1)+diag(b)+diag(c,1);
D = tridiag;
Nframes=5;
marca=floor(tsteps/(Nframes-1));
u_save = zeros(2*M+1,Nframes);
%
t_save = zeros(1,Nframes);
%f_save = zeros(2*M+1,Nframes);
%f_save =
%Le ponemos las condiciones de frontera
u_save(1,:)=zeros(1,Nframes);
u_save(2*M+1,:)=zeros(1,Nframes);
%guardamos la posicion de partida
u_save(2:2*M,1)=u;
t_save(1)=0;
%Bucle principal
I=eye(2*M-1);
%
A=(I+dt*D);
%
% método 1 es para ver la supeficie y el método 2 es para ver la solución
% para algunos tiempos
%metodo=METODO;
%switch metodo
% case 1
%superficie
% u_save(1,:)=zeros(1,tsteps);
% u_save(2*N+1,:)=zeros(1,tsteps);
% for n=1:tsteps
% B=fi(2:2*N,n);
% u=A*u+dt*B;
% u_save(2:2*N,n)=u;
%
% end
%
% %
% figure(2)
% mesh(u_save)
%
% case 2
u = A*u+dt*fi(2:2*M,1);
for n=1:tsteps
B=fi(2:2*M,n);
u=A*u+dt*B;
% u_save(2:2*N,n)=u;
%end
% for n=1:tsteps
% u=A*u+b;
%Guardamos los valores de u para algunos tiempos
if mod(n,marca)==0
indice=1+n/marca;
u_save(2:2*M,indice)=u;
t_save(indice)=tiempos(n);
end
end
figure(2)
plot(x,u_save(:,3))
%figure(3)
%plot(x,u_save)
%end
1 commentaire
Walter Roberson
le 14 Avr 2018
What assistance are you requesting?
Réponses (2)
Olivar Luis Eduardo
le 16 Avr 2018
0 votes
1 commentaire
Walter Roberson
le 16 Avr 2018
Are you getting an error message? What difference are you seeing between what you expect and what you get? What are the differential equations, in symbolic form (such as an image of the equations) ?
Olivar Luis Eduardo
le 16 Avr 2018
0 votes
1 commentaire
Walter Roberson
le 16 Avr 2018
I do not understand about Fallanda ? I find reference to a breed of horses, but I do not see anything that might relate to this question??
Catégories
En savoir plus sur Symbolic Solvers dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!