• API
• FAQ
• Tools
• Archive
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.
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)
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.

Top