RK4 and Modeling Radiovirotherapy
Fall 2021
1. Instructions
This is the semester project and is due by midnight on the day of the final exam. I encourage
you all to work together on it, but do not simply copy one another. You must turn in all your
code and plots via Canvas.
2. RK4 Algorithm
Assume that the solution u(t) is approximated by y(t). That is, yj ≈ u(tj). To approximate the
solution of the nth-order system of first-order IVPs given by
du1
dt
= f1(t, u1, u2, . . . , un)
du2
dt
= f2(t, u1, u2, . . . , un)
...
dun
dt
= fn(t, u1, u2, . . . , un)
for t ∈ [α, β] and initial data
u1(α) = a1
...
un(α) = an,
for (N + 1) uniformly spaced nodes in [α, β], use the following algorithm.
(1) Define the step size by h = β−α
N
and set t = α.
(2) For j = 1, . . . n, define yj = aj.
(3) For current t, compute the following for each yj with j = 1, . . . , n
k1,j = hfj (t, y1, . . . , yn)
k2,j = hfj
(
t+
h
2
, y1 +
1
2
k1,1, y2 +
1
2
k1,2, . . . , yn +
1
2
k1,n
)
k3,j = hfj
(
t+
h
2
, y1 +
1
2
k2,1, y2 +
1
2
k2,2, . . . , yn +
1
2
k2,n
)
k4,j = hfj (t+ h, y1 + k3,1, y2 + k3,2, . . . , yn + k3,n)
Update yj+1:
yj+1 = yj +
1
6
(k1,j + 2k2,j + 2k3,j + k4,j)
(4) Update t.
Create a solver and run file that will successfully implement the RK4 method for systems of
differential equations for n = 2 and n = 3. Your solver should do the following:
(1) The single solver should have an input option to run either RK4 with n = 2 or n = 3. I
recommend you code them individually first and then combine them into a single solver.
(2) Plot the approximate solution for each yj and the exact solution u(tj) on the same set
of axes.
(3) Show a table in the command window that records the absolute error at each t-value for
each uj when the exact solution is known. Otherwise, it should just show a table in the
command window that shows the yj value at each t-value.
3. Part 1: Model Problems
(1) Run your code on the following 2 test codes. For the n = 2 case, numerically approximate
the solution to the differential equation given below.
u′1 = −4u1 + 3u2 + 6
u′2 = −2.4u1 + 1.6u2 + 3.6
with u1(0) = 0 and u2(0) = 0. The exact solutions are given by
u1(t) = −3.375e−2t + 1.875e−0.4t + 1.5
u2(t) = −2.25e−2t + 2.25e−0.4t
Use N = 5. Check that your results agree with the table below.
tj y1,j y2,j |u1(tj)− y1,j| |u2(tj)− y2,j|
0.000000 0.000000 0.000000 0.000000 0.000000
1.00000E-01 5.38260E-01 3.19630E-01 8.70680E-06 5.80370E-06
2.00000E-01 9.68500E-01 5.68780E-01 1.42570E-05 9.50280E-06
3.00000E-01 1.31070E+00 7.60730E-01 1.75080E-05 1.16700E-05
4.00000E-01 1.58130E+00 9.06320E-01 1.91110E-05 1.27380E-05
5.00000E-01 1.79350E+00 1.01440E+00 1.95580E-05 1.30350E-05
For the n = 3 case, numerically approximate the solution to the differential equation
given below.
u′1 = u2
u′2 = −u1 − 2et + 1
u′3 = −u1 − et + 1
with u1(0) = 1, u2(0) = 0, and u3(0) = 1 for 0 ≤ t ≤ 2 and h = 0.5. The actual solutions
are given by
u1(t) = cos t+ sin t− et + 1
u2(t) = − sin t+ cos t− et
u3(t) = − sin t+ cos t.
Check that your results agree with the table below.
tj y1,j y2,j y3,j |u1(tj)− y1,j| |u2(tj)− y2,j| |u3(tj)− y2,j|
0 1 0 1 0 0 0
5.00000E-01 7.07870E-01 -1.24990E+00 3.98850E-01 4.16070E-04 6.77620E-04 6.91590E-04
1.00000E+00 -3.36920E-01 -3.01760E+00 -2.99320E-01 4.08990E-04 1.80870E-03 1.84570E-03
1.50000E+00 -2.41330E+00 -5.40520E+00 -9.23470E-01 1.29540E-04 3.21410E-03 3.28910E-03
2.00000E+00 -5.89480E+00 -8.70970E+00 -1.32050E+00 1.11540E-03 4.79500E-03 4.93260E-03
4. Part 2: Application to Cancer Radiovirotherapy*
This portion of the project is subject to variability for each student. If you would like to apply
your RK4 code to some other type of application, contact your instructor prior to 1 month
before the due date (midnight of the day of the final). Otherwise, you will complete (1) through
(3) below.
(1) Read the paper titled Mathematical modeling of cancer radiovirotherapy by Dingli et al.
(2) Write 3-4 paragraph summary of the work. You do not need to understand everything
in the paper, but give the “big picture” idea.
(3) Use your solver to recreate figures 5 and 6. You will need to sift through the paper to
find the parameters used in their work. Your plots should include appropriate titles and
legends.

欢迎咨询51作业君