kill(all)$

appendix : electrical  , solve two differential equation , numerical & symbolic,version 1.0,part2
peter.vlasschaert@gmail.com, 11/01/2020

part1 :symbolic       ( second differential  equatio ,see 'below'        )., use p2c
==============
part 2 : numerical    (system of two first order differential equations: newton ,rlc -circuit).
==============

load packages
---------------------

load(lrats)$

solve problem 'part 1' ,numerical
=========================

1e)  second order differential equations
-------------------------------------------------------

declare(k,constant);

p1:'diff(y,x,2)=-k*x;

2e)  1e by two first order equations.
---------------------------------------------------

p2a: 'diff(y,x) = z;
p2b: 'diff(z,x)  = -k*x;

p2a1: rhs(p2a);
p2b1:rhs(p2b);

repeat :part 1

p2c1:subst([k=1],p1);
p2c2:ode2(p2c1,y,x);
p2c3:ic2(p2c2, x=0, y=0, 'diff(y,x)=1);
p2c4:rhs(p2c3);

3e) create a function rhs 2e)
---------------------------------------

k:1$

F(x,y,z):=''p2a1;
G(x,y,z,k) := ''p2b1;

4e) system  of ode equations 'p2a,p2b',rk4
-----------------------------------------------------------

n:20;
x0:0;
y0:0;
z0:1;
h:0.25;

 vx: makelist(0, n + 1);   /* reserved array */
 vy: makelist(0, n + 1);  
 vz: makelist(0, n + 1);  
   x: x0; 
   y: y0;
   z:z0;
   vx[1]: x0;
   vy[1]: y0;
   vz[1]:z0;
   for i from 1 thru n do (
      f1: ev(h * F(x, y,z),float),
       g1: ev(h * G(x, y,z,k),float),
      f2: ev(float(h * F(x + h / 2, y + f1 / 2,z+f1/2))),
       g2: ev(float(h * G(x + h / 2, y + g1 / 2,z+g1/2,k))),
      f3: ev(float(h * F(x + h / 2, y + f2 / 2,z+f2/2))),
        g3: ev(float(h *G(x + h / 2, y +g2 / 2,z+g2/2,k))),
      f4: ev(float(h * F(x + h, y + f3,z+f3))),
        g4: ev(float(h *G(x + h, y + g3,z+g3,k))),
      y: y + (f1 + 2 * f2 + 2 * f3 + f4) / 6,
      vy[i + 1]: y,
      z: z + (g1 + 2 * g2 + 2 * g3 + g4) / 6,
      vz[i + 1]: z,
      x:x+h ,         
      vx[i + 1]: x)$

vx;
vy;
vz;

plot : result 'ode'
============

1e) plot function vx versus vy
======================

x;
y;

kill(x,y);

x;
y;

reason : plot option to kill  values x,y

v:[]$
for i from 1 thru n+1 do(
v:endcons([vx[i],vy[i]],v))$
v;
wxplot2d([discrete,v])$
p2c:length(v);

2e) exact solution: part 1: short version above,'p2c4'
========================================

p3:p2c4;

menu : Plot :Plot2d

wxplot2d([%], [x,0,5.05])$

x;
y;

 3e) both exact &  numerical approximation : same plot
==========================================

p4:[[discrete,v],p3]$

wxplot2d(%,[x,0,5.05])$

4e)discrete values :points, exact solution :lines
====================================

p4a:[style,[points,2,2.6],[lines,1.5]]$
p4b:[x,0,5]$
p4c1:[xlabel,"x values"]$
p4c2:[ylabel,"y"]$
p4c3:[legend,"discrete","exact"]$
wxplot2d(p4,p4b,p4a,p4c1,p4c2,p4c3)$

                    ================== <  plot  > ============================

kill(x,y,z);

[x,y,z];

how rk from ' Maxima' to solve
=======================

1e)  load packages
===============

load(diffeq)$

2e) use rk from 1e)
===============
model                      :         =>  'diff(y,x,2)=-k*x
(I)nitial (C)ondition :       1e) x = 0  => y =0
                                         2e) x = 0  => diff(y,x) = 1
see: also exact solution above.

a) model:
-------------

p5a:p2a;
p5b:'diff(z,x) =-x;

see:google

http://maxima.sourceforge.net/docs/manual/maxima_22.html#Item_003a-rk

==========================================================
 dx/dt = 4-x^2-4*y^2    
 dy/dt = y^2-x^2+1
(%i1) sol: rk([4-x^2-4*y^2,y^2-x^2+1],[x,y],[-1.25,0.75],[t,0,4,0.02])$
(%i2) plot2d ([discrete,makelist([p[1],p[3]],p,sol)], [xlabel,"t"],[ylabel,"y"])$
==========================================================

b)  solution & plot 'rk'
-----------------------------

sol:rk([rhs(p5a),rhs(p5b)],[y,z],[0,1],[x,0,5,h])$
p6:makelist([p[1],p[2]],p,sol);
p6a1:length(p6);
wxplot2d ([discrete,p6], [xlabel,"x values"],[ylabel,"y"])$

combination : three plots
===================

plot : 'rk4,exact,rk ' versus x
=====================

p6:[[discrete,v],p3,[discrete,p6]]$
p6a:[style,[points,2,2.6],[lines,1],[points,2,1.6]]$
p6b:p4b$
p6c:p4c1$
p6d:p4c2$
p6e:[legend,"rk4","exact","rk"]$
wxplot2d(p6,p6b,p6a,p6c,p6d,p6e)$

[x,y,z];

reason : x,y,z ,var

construct definition : rk4 , use ' block'
============================

kill(f1,f2,f3,f4)$
kill(g1,g2,g3,g4)$

r2k4(F,G, a1,  a2, a3,h, n) := block(
      x:a1,
      y:a2,
      z:a3,
      f1: ev(h * F(x, y,z),float),
       g1: ev(h * G(x, y,z,k),float),
      f2: ev(float(h * F(x + h / 2, y + f1 / 2,z+f1/2))),
       g2: ev(float(h * G(x + h / 2, y + g1 / 2,z+g1/2,k))),
      f3: ev(float(h * F(x + h / 2, y + f2 / 2,z+f2/2))),
        g3: ev(float(h *G(x + h / 2, y +g2 / 2,z+g2/2,k))),
      f4: ev(float(h * F(x + h, y + f3,z+f3))),
        g4: ev(float(h *G(x + h, y + g3,z+g3,k))),
      y: y + (f1 + 2 * f2 + 2 * f3 + f4) / 6,
      vyy[i + 1]: y,
      z: z + (g1 + 2 * g2 + 2 * g3 + g4) / 6,
      vzz[i + 1]: z,
      x:x+h ,         
      vxx[i + 1]: x)$

   vxx: makelist(0, n + 1)$
   vyy: makelist(0, n + 1)$
   vzz: makelist(0, n + 1)$
   vxx[1]: x0;
   vyy[1]: y0;
   vzz[1]: z0;
   for i from 1 thru n  do
                                     (
      r2k4(F,G, vxx[i],  vyy[i] ,vzz[i] ,h, n) 
                                     )$

vxx;
vyy;
vzz;

plot :use definition
==============

vv:[]$
for i from 1 thru n+1 do(
vv:endcons([vxx[i],vyy[i]],vv))$
vv$

p7:[[discrete,vv],p3]$

p7a:wxplot2d(p7,p6b,[legend,"discrete","exact"],[xlabel," x-values"],[ylabel,"y"])$