Help me please with this code:

Question

Olivar Luis Eduardo 2018-4-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/394949-help-me-please-with-this-code

评论： Walter Roberson 2018-4-16

%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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2018-4-14

What assistance are you requesting?

http://www.mathworks.com/matlabcentral/answers/6200-tutorial-how-to-ask-a-question-on-answers-and-get-a-fast-answer

请先登录，再进行评论。

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

Answer 1

Olivar Luis Eduardo 2018-4-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/394949-help-me-please-with-this-code#answer_315470

I can not reconstruct the solution with this code, what am I failing?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2018-4-16

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) ?

请先登录，再进行评论。

Answer 2

Olivar Luis Eduardo 2018-4-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/394949-help-me-please-with-this-code#answer_315471

This code is used to find the solution of the parabolic equation in partial differential equations, using finite differences. Where am I from Fallanda?