RK4 fungerer

finitedifferences.py

@@ -4,11 +4,24 @@ def forwardEuler(m):
     for i in range(0,m):
         if(i==m-1):
             A[i,i]=-1
-            A[i,i-1]=0
-            A[i,0]=1
+            A[i,i-1]=1
+            A[i,0]=0
         else:
             A[i,i]=-1
-            A[i,i-1]=0
-            A[i,i+1]=1
+            A[i,i-1]=1
+            A[i,i+1]=0
 
-    return A
+    return A
+
+def centralDifference(m):
+    A = np.zeros([m,m])
+    for i in range(0,m):
+        if(i==m-1):
+            A[i,i]=0
+            A[i,i-1]=-0.5
+            A[i,0]=0.5
+        else:
+            A[i,i]=0
+            A[i,i-1]=-0.5
+            A[i,i+1]=0.5
+    return A

test.py

@@ -17,19 +17,31 @@ def solver(A,f,u,h,k):
 
     return u
 
+def RK4(A,f,u,h,k):
+    u[:,0]=f[:]
+    for i in range(1,u.shape[1]):
+        k1 = k/h*A @ u[:,i-1]
+        k2 = k/h*A @ (u[:,i-1] + 0.5*k1)
+        k3 = k/h*A @ (u[:,i-1] + 0.5*k2)
+        k4 = k/h*A @ (u[:,i-1] + k3)
+        u[:,i] = u[:,i-1] + 1/6*(k1 + 2*k2 + 2*k3 + k4)
+    return u
+
 
 
-T=1
-m=300
-n = 3590
+T=2
+m=15
+n = 50
 u = np.zeros([m,n])
 h=1/m
 k = T/n
-A = finitedifferences.forwardEuler(m)
+A = finitedifferences.centralDifference(m)
+
+
 x = np.linspace(0,1,m)
 f = np.sin(2*3.1415*x)
-u[:,0] = f[:]
-u = solver(A,f,u,h,k)
+print(u.shape[1])
+u = RK4(A,f,u,h,k)
 X,Y=np.meshgrid(np.linspace(0,1,m),np.linspace(0,T,n))
 print(X.shape)
 print(Y.shape)
@@ -40,4 +52,7 @@ def solver(A,f,u,h,k):
 ax = plt.axes(projection='3d')
 
 ax.plot_surface(X, Y, u.T,cmap='viridis', edgecolor='none')
+
+ax.set_xlabel('X axis')
+ax.set_ylabel('t axis')
 plt.show()