SHARE
TWEET

Untitled

a guest Dec 6th, 2019 80 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. clear all
  2. mu=0.01; tau=100; L=1; M=100; N=2000; x=linspace(0,L,M)';
  3. v=sqrt(tau/mu); dx=x(2)-x(1); dt=dx/v; p=(v*dt/dx)^2;
  4. % s=-2*x.^2 + 0.8*x + 0.92;
  5. % g=transpose(sin((1:M)*2*3.14));  
  6. % r=-0.28*ones(length(1),N);
  7. % h=0.8 * ones(length(1),N);
  8. s=(x+1/9).^(1/2)-0.5;
  9. g=1*exp(-(x-1/2).^2)-1*exp(-0.25);
  10. l=-1/6*ones(1,N);
  11. k=3*sqrt(10)/20*ones(1,N);
  12.  
  13. b = 0.5;
  14. beta = b/(2*mu);
  15. q = 1+beta*dt;
  16. u = 1-beta*dt;
  17.  
  18. f(:,1)=s;
  19. f(1,2)=l(2);
  20. f(2:M-1,2)=p/2*(f(3:M,1)-2*f(2:M-1,1)+ f(1:M-2,1))+f(2:M-1,1)+u*dt*g(2:M-1);
  21. f(M,2)=p*(dx*k(1)+f(M-1,1)-f(M,1))+f(M,1)-u*dt*g(M);
  22. for n=2:N-1
  23. f(1,n+1)=l(n+1);
  24. f(2:M-1,n+1)=p/q*(f(3:M,n)-2*f(2:M-1,n)+f(1:M-2,n))+2/q*f(2:M-1,n)-u/q*f(2:M-1,n-1);
  25. f(M,n+1)=2*p/q*(dx*k(n)+f(M-1,n)-f(M,n))+2/q*f(M,n)-u/q*f(M,n-1);
  26. clf;
  27. plot(1:M,f(:,n));
  28. ylim([-1,1]);
  29. drawnow;
  30. end
  31.  
  32. %zadanie przykladowe
  33. % f(:,1)=s;
  34. % f(1,2)=p*(f(2,1)- f(1,1)-dx*h(1)) + f(1,1) +u*dt*g(1);
  35. % f(2:M-1,2)=p/2*(f(3:M,1)-2*f(2:M-1,1)+ f(1:M-2,1))+f(2:M-1,1)+u*dt*g(2:M-1);
  36. % f(M,2)=r(2);
  37. % for n=2:N-1
  38. % f(1,n+1)=2*p/q*(f(2,n)-f(1,n)- dx*h(n) )+2/q*f(1,n)-u/q*f(1,n-1);
  39. % f(2:M-1,n+1)=p/q*(f(3:M,n)-2*f(2:M-1,n)+f(1:M-2,n))+2/q*f(2:M-1,n)-u/q*f(2:M-1,n-1);
  40. % f(M,n+1)=r(n+1);
  41. % end
  42.  
  43.  
  44. % f(:,1)=s;
  45. % f(1,2)=p*(f(2,1)- f(1,1)-dx*h(1)) + f(1,1) +dt*g(1);
  46. % f(2:M-1,2)=p/2*(f(3:M,1)-2*f(2:M-1,1)+ f(1:M-2,1))+f(2:M-1,1)+dt*g(2:M-1);
  47. % f(M,2)=r(2);
  48. % for n=2:N-1
  49. %     f(1,n+1)=2*p*(f(2,n)-f(1,n)- dx*h(n) )+2*f(1,n)-f(1,n-1);
  50. %     f(2:M-1,n+1)=p*(f(3:M,n)-2*f(2:M-1,n)+f(1:M-2,n))+2*f(2:M-1,n)-f(2:M-1,n-1);
  51. %     f(M,n+1)=r(n+1);
  52. %     clf;
  53. %     plot(1:M,f(:,n));
  54. %     ylim([-4,4]);
  55. %     drawnow;
  56. % end
  57.  
  58. surf(f)
  59. shading interp
  60.  
  61. %equilibrium
  62. % plot(x, f(:,N),'LineWidth', 2);
  63. %  
  64. % hold on
  65. % yE = 0.8.*x - 1.08;
  66. % plot(x, yE,'r','LineWidth', 2);
  67. % grid on
  68. %  
  69. % legend('last iteration', 'equilibrium');
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top