copy the past then run the matlab program with varying dx and dt for multiple ex

Question

copy the past then run the matlab program with varying

dx and dt for multiple experiments then explain the

results.

clear all

clc

close all

   c=20.;a=10.;pi=3.14;

   dx=20000.;dt=100.;

% c m: number of time steps; n: number of grid points

   m=5000;n=50

   u1(1:m,1:n)=0;

   ek1(1:m)=0;

% c initial condition

   l=10.*dx;

   kk=2*pi/l;

   for j=1:n

u1(1,j)=c+a*sin(kk*dx*j);

   end

% c boundary Condition (static)

   for i=1:m

u1(i,1)=u1(1,1);

u1(i,n)=u1(1,n);

   end

% c Caculate u1: leapfrog scheme

% c first-step: upstream scheme

for j=2:n-1

u1(2,j)=u1(1,j)-c*dt*(u1(1,j)-u1(1,j-1))/dx;

end

% % % c time integration: leapfrog scheme

   for i=3:m

   for j=2:n-1

   u1(i,j)=u1(i-2,j)-c*dt*(u1(i-1,j+1)-u1(i-1,j-1))/dx;

   end

   end

% % c time integration: upstream scheme

% for i=3:m

% for j=2:n-1

% u1(i,j)=u1(i-1,j)-c*dt*(u1(i-1,j+1)-u1(i-1,j-1))/dx;

% end

% end

figure(1)

   plot (1:50,u1(m,:),'+');

   hold on

   plot(1:50,u1(m,:),'b');

   xlabel ('grid points');

   ylabel ('u');

% c calulate total krinetic energy for every time step

   for i=1:m

   ek1(i)=0.;

   end

   for i=1:m

   for j=1:n

ek1(i)=ek1(i)+0.5*u1(i,j)^2;

   end

   end

   mm=1:m

   figure(2)

   plot (mm(1:5000),ek1(1:5000),'+');

   xlabel ('Time steps');

   ylabel ('Energy')

% c output

% c print KE for every 50 time steps

   fid=fopen('nwpc-200.dat','w');

for i=1:50:m

if mod((i-1)/50,5)==4

fprintf(fid,' %10.2E ',ek1(i));

else

fprintf(fid,' %10.2E ',ek1(i));

end

end

   fclose(fid)

Goals: 1) Practice numerical methods to solve simple PDEs with finite difference equations; 2) Understand the linear computational stability Problem # 1 Integrate the linear wave equation using values typical of large-scale models. You can write your own MATLAB program or modify a sample program I offered to you (see the attachment). au Boundary conditions: periodic Initial conditions: ux,0)c+Asin(kx) c = 20ms-1, A-10 ns" ,Ar = 200km, k = 2 / L withL = 10Ax - (a) Use leapfrog scheme. Choose two time steps: one satisfies the CFL condition and the other violates it. 1) How long does it take to "blow up"; 2) Calculate the total kinetic energy of the solution u (0.5 u') at each time step; write a simple MATLAB (or IDL) program and draw a curve to show the variation of total kinetic energy as function of the number of time steps. (b) Compare with the exact solution, computer the root-mean-square (RMS) error R and the relative error RE. Then, Repeat with A-25 m/s Repeat with L = 4Dx Prepare a table to summarize R and RE. (c) Modify the programs; use the upstream scheme (instead of leapfrog scheme) and repeat (a) Lab Report: Based on your results to solve above problem, write a report (minimum 1 to 2 pages of text for each problem; 12pt Times New Roman; doubled space) to discuss the linear computational stability problem. You can attach figures to support your conclusion.

Explanation / Answer

clear all
clc
close all
c=20.;a=10.;pi=3.14;
dx=20000.;dt=100.;
% c m: number of time steps; n: number of grid points
m=5000;n=50

u1(1:m,1:n)=0;
ek1(1:m)=0;

% c initial condition
l=10.*dx;
kk=2*pi/l;
for j=1:n
u1(1,j)=c+a*sin(kk*dx*j);
end

% c boundary Condition (static)
for i=1:m
u1(i,1)=u1(1,1);
u1(i,n)=u1(1,n);
end

% c Caculate u1: leapfrog scheme
% c first-step: upstream scheme
for j=2:n-1
u1(2,j)=u1(1,j)-c*dt*(u1(1,j)-u1(1,j-1))/dx;
end
% % % c time integration: leapfrog scheme
for i=3:m
for j=2:n-1
u1(i,j)=u1(i-2,j)-c*dt*(u1(i-1,j+1)-u1(i-1,j-1))/dx;
end
end
% % c time integration: upstream scheme
% for i=3:m
% for j=2:n-1
% u1(i,j)=u1(i-1,j)-c*dt*(u1(i-1,j+1)-u1(i-1,j-1))/dx;
% end
% end

figure(1)
plot (1:50,u1(m,:),'+');
hold on
plot(1:50,u1(m,:),'b');
xlabel ('grid points');
ylabel ('u');

% c calulate total krinetic energy for every time step
for i=1:m
ek1(i)=0.;
end

for i=1:m
for j=1:n
ek1(i)=ek1(i)+0.5*u1(i,j)^2;
end
end

mm=1:m
figure(2)
plot (mm(1:5000),ek1(1:5000),'+');
xlabel ('Time steps');
ylabel ('Energy')

% c output
% c print KE for every 50 time steps
fid=fopen('nwpc-200.dat','w');
for i=1:50:m
if mod((i-1)/50,5)==4
fprintf(fid,' %10.2E ',ek1(i));
else
fprintf(fid,' %10.2E ',ek1(i));
end
end

fclose(fid)

Related Questions

Navigate

• Browse (All)
• Browse C
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.