From 1d8d90c26c7780b91f5f3a681fbb55c51642a96a Mon Sep 17 00:00:00 2001
From: SorBern <64829934+SorBern@users.noreply.github.com>
Date: Sun, 10 May 2020 01:29:48 +0200
Subject: [PATCH] RK4 and Euler for different f()
---
fluids4.py | 143 +++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 143 insertions(+)
create mode 100644 fluids4.py
diff --git a/fluids4.py b/fluids4.py
new file mode 100644
index 0000000..50c8e47
--- /dev/null
+++ b/fluids4.py
@@ -0,0 +1,143 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun May 10 01:11:03 2020
+
+@author: Marina
+"""
+
+from numpy import *
+
+def RK4(Xstart,Ustart,Xend,h,f):
+ imax = int((Xend-Xstart)/h)
+ X = Xstart + arange(imax+1)*h
+ U = zeros((size(Ustart), imax+1))
+ U[:, 0] = asarray(Ustart)
+ for i in range(1,imax+1):
+ K1 = asarray(f(U[:, i-1]))
+ K2 = asarray(f(U[:, i-1] + (h/2) * K1))
+ K3 = asarray(f(U[:, i-1] + (h/2)* K2))
+ K4 = asarray(f(U[:, i-1] + h * K3))
+
+ U[:, i] = U[:, i-1] + (h/6) * (K1 + 2 * K2 + 2 * K3 + K4)
+
+ U = transpose(U)
+ return X, U
+
+
+def Euler(Xstart,Ustart,Xend,h,f):
+ N = 5
+ X = arange(Xstart,Xend,h)
+ n = size(X)
+ U = zeros((N, n))
+ U[:, 0] = transpose(Ustart)
+ for j in range (1, n):
+ U[:,j] = U[:, j-1] + h*f(U[:,j-1])
+
+ U = transpose(U)
+ return X,U
+
+# We define:
+# f = u1; f' = v1; f''= w1
+# \u03B8 = u2; \u03B8'= v2
+# u = [u1, v1, w1, u2, v2]
+
+def f(u):
+ du1dt = u[2]
+ du2dt = u[3]
+ dv1dt = u[4]
+ dv2dt = -3*0.01*u[0]*u[3]
+ dw1dt = -3*u[0]*u[4] + 2*u[2]**2 - u[1]
+ return array([du1dt,dv1dt,dw1dt,du2dt,dv2dt])
+
+
+# Shoot method
+
+def shoot(a,b,h,integrator):
+ X,U = integrator(0,[0,1,0,a,b],30,h,f)
+ print(-U)
+ return array([-U[-1,3], -U[-1,4]])
+
+# We solve the problem in the same way we solved the Blasius equation,
+# aplaying the bisection method.
+# We are supposing that we want to get the an error equal to the tolerance in both cases,
+# to calculate alfa1 and alfa2, and that we can choose different intervals.
+
+def blasius(delta1,delta2,nmax,tol,h,integrator):
+ delta0 = array([delta1, delta2])
+ a = zeros(2)
+ b = zeros(2)
+ a[0] = delta0[0,0]
+ a[1] = delta0[1,0]
+ b[0] = delta0[0,1]
+ b[1] = delta0[1,1]
+ x = zeros(2)
+ delta = zeros(2)
+
+ print('a1,a2',a[0],a[1])
+ fa = shoot(a[0],a[1],h,integrator)
+ print('fa', fa)
+ print('b1,b2',b[0],b[1])
+ fb = shoot(b[0],b[1],h,integrator)
+ print('fb', fb)
+ n = 1
+ delta[0] = (b[0]-a[0])/2
+ delta[1] = (b[1]-a[1])/2
+ for i in range (0,2):
+ if (fa[i]*fb[i] > 0) :
+ return 0, 0,'Bad initial interval :-( for i = %d' % (i)
+
+ while (abs(delta[0]) >= tol and abs(delta[1]) >= tol and n < nmax) :
+ n = n + 1
+ for i in range (0,2):
+ delta[i] = (b[i]-a[i])/2
+ x[i] = a[i] + delta[i]
+ print('x%d = %13.7e'%(i,x[i]))
+
+
+ fx = shoot(x[0],x[1],h,integrator)
+ print('fx', fx)
+ print(" x1 = %14.7e (Estimated error %13.7e at iteration %d)" % (x[0],abs(delta[0]),n))
+ print(" x2 = %14.7e (Estimated error %13.7e at iteration %d)" % (x[1],abs(delta[1]),n))
+
+ for i in range (0,2):
+ if (fx[i]*fa[i] > 0) :
+ a[i] = x[i]
+ fa[i] = fx[i]
+ else:
+ b[i] = x[i]
+ fb[i] = fx[i]
+
+ if (n == nmax) :
+ return x[0], x[1],'Increase nmax : more iterations are needed :-('
+ return x[0], x[1],'Convergence observed :-)'
+
+# These are the initial values we have choosen
+
+Xend = 30
+h = 0.5
+Pr = 0.01
+delta1 = [-1,0]
+delta2 = [1,1.5]
+nmax = 20
+tol = 1e-4
+
+a, b, message = blasius(delta1,delta2,nmax,tol,h,RK4)
+X,U = RK4(0,[0,1,0,a,b],Xend,h,f)
+print(U,shape(U))
+
+print('For RK4:')
+print(" === Requested f''(0) = %.4f === %s" % (a,message))
+print(" === Requested \u03B8'(0) = %.4f === %s" % (b,message))
+print(" === Obtained final value for f'(0) = %13.7e " % U[-1,3])
+print(" === Obtained final value for \u03B8(0) = %13.7e " % U[-1,4])
+
+from matplotlib import pyplot as plt
+import matplotlib
+matplotlib.rcParams['toolbar'] = 'None'
+plt.rcParams['figure.facecolor'] = 'lavender'
+plt.rcParams['axes.facecolor'] = 'lavender'
+
+fig = plt.figure("Fluids")
+plt.plot(X,U[:,3],'-b',X,U[:,4],'-r')
+plt.text(4,2e48,"f'(x)",color='blue',fontsize=12)
+plt.text(0.5,1e48,"\u03B8 (x)",color='red',fontsize=12)
\ No newline at end of file
--
GitLab