Trinh @ Bath

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Next revision
Previous revision
ara-prob5 [2021/03/23 08:05]
trinh created
ara-prob5 [2021/03/26 11:29] (current)
trinh
Line 1: Line 1:
 +====== Monotonic shock solutions ======
 +
 +This is the caller code. You can call it 'KS_monotonic_caller.m'. Copy the below script and paste the code in and save it in an appropriate directory.
 +
 <code> <code>
 % KS_MONOTONIC_CALLER will solve the K-S equation for the case of the % KS_MONOTONIC_CALLER will solve the K-S equation for the case of the
Line 36: Line 40:
 drawnow drawnow
 </code> </code>
 +
 +The caller code will require 'KSode.m'. Save it and run it from the same directory. Make sure you name it with the correct upper and lower-case characters.
  
 <code> <code>
Line 50: Line 56:
 Yp = [up; upp; uppp]; Yp = [up; upp; uppp];
 </code> </code>
 +
 +====== Oscillatory shock solutions ======
 +
 +Again, you can copy and paste the below script into a Matlab file and call it, e.g. 'KS_osc_caller.m'. It will make use of the 'KSode.m' code above.
 +
 +<code>
 +% KS_OSC_CALLER will solve the K-S equation for the case of the
 +% oscillatory shock conditions
 +% -------------------------------------------------------------------
 +%
 +%   Written 19 Mar 2021 for the INI Spring School
 +%           EXPONENTIAL ASYMPTOTICS FOR PHYSICAL APPLICATIONS
 +
 +z0 = 1.24;
 +phi = 2.60586555;
 +% phi = 2.59;
 +% phi = 2.5;
 +
 +ep = 0.05;
 +zmin = -6;
 +zmax = 12;
 +
 +gam = sqrt(1/ep^2 - 1);
 +bc = @(z) -1 + exp(-(z-z0))*sin(gam*(z-z0)+phi);
 +
 +bcp = @(z) -exp(-(z-z0))*sin(gam*(z-z0)+phi) + ...
 +    gam*exp(-(z-z0))*cos(gam*(z-z0)+phi);
 +
 +bcpp = @(z) exp(-(z-z0))*sin(gam*(z-z0)+phi) ...
 +    - 2*gam*exp(-(z-z0))*cos(gam*(z-z0)+phi) ...
 +    -gam^2*exp(-(z-z0))*sin(gam*(z-z0)+phi);
 +ic = [bc(zmax); ...
 +    bcp(zmax); ...
 +    bcpp(zmax)];
 +
 +
 +fwd = @(t,Y) KSode(t,Y,ep);
 +
 +mytol = 1e-10;
 +options = odeset('RelTol', mytol, 'AbsTol', mytol);
 +sol = ode113(fwd, [zmax, zmin], ic, options);
 +
 +figure(1)
 +hold all
 +plot(sol.x, sol.y(1,:));
 +ylim([-5,5]);
 +drawnow
 +</code>
 +