DjangoBB LoFi version

Полная версия: Вопрос по Рунге-Кутта 4(5) (Метод Рунге-Кутта-Мерсона)

Начало » Центр помощи » Вопрос по Рунге-Кутта 4(5) (Метод Рунге-Кутта-Мерсона)

Wahlberg

Дек. 6, 2012 12:52:34

Собственно есть вот такое задание

Задание 9

Построить дифференциальные уравнения плоского движения твёрдого тела под действием постоянной силы.

Построить функцию численного интегрирования дифференциальных уравнений методом Рунге-Кутта 4(5) (Метод Рунге-Кутта-Мерсона). Проинтегрировать уравнения движения твёрдого тела.

Решая данную задачу получаем СДУ(система дифференциальных уравнений):

(V_x)'=F*cos(phi)/m
(V_y)'=F*sin(phi)/m
Omega'=F*a/J
x'=V_x
y'=V_y
Phi'=omega

Собственно вот мой код.

import math
import pylab
from matplotlib import mlab
from scipy.linalg import solve
from scipy import array
import matplotlib.pyplot as plt 
from pylab import *

m=10
F=100
J=10
a=20
count=F*a/J
Nu=F/m

def f(t,var):
    return t*count

def rkm45(h,t,t0,tmax,var):
    i=0
    item=zeros((tmax-t)/h)
    while t0<tmax:
        t0=t0+h
        t=t+h
        k1=f(t,var)
        k2=f(t+h/3,var+h/3*k1)
        k3=f(t+h/3,var+h/6*k1+h/6*k2)
        k4=f(t+h/2,var+h/8*k1+3*h/8*k2) 
        k5=f(t+h,var+h/2*k1-3*h/2*k3+2*h*k4)    
        item[i]=var
        i=i+1        
        var=var+h*(k1+4*k4+k5)/6
    print item
    return item
    
rkm45(0.05,0,0,2,0)

Собственно вопрос, как задать функцию от функции, если то что выше дай бог правильно написано.

ziro

Дек. 7, 2012 07:34:00

В питоне функции похожи на обычные переменные. Вы можете смело ее передавать как параметр другой функции. Ну что-то типа такого:

def rkm45(f, h,t,t0,tmax,var):
    i=0
    item=zeros((tmax-t)/h)
    while t0<tmax:
        t0=t0+h
        t=t+h
        k1=f(t,var)
        k2=f(t+h/3,var+h/3*k1)
        k3=f(t+h/3,var+h/6*k1+h/6*k2)
        k4=f(t+h/2,var+h/8*k1+3*h/8*k2) 
        k5=f(t+h,var+h/2*k1-3*h/2*k3+2*h*k4)    
        item[i]=var
        i=i+1        
        var=var+h*(k1+4*k4+k5)/6
    print item
    return item
    
rkm45(f, 0.05,0,0,2,0)

Wahlberg

Дек. 11, 2012 01:11:59

Сделал, вдруг кому понадобится

))

# -*- coding: utf-8 -*-
"""
Created on Sun Dec 09 18:28:40 2012

@author: ваха
"""
    
from numpy import *
import pylab
from matplotlib import mlab
from pylab import *

a=2
b=3
m=1
F=13
t0=0
tmax=10
epsilon=0.01
y0=[0,0,math.pi/3,0,0,0]
h=0.1

def rkm45(f,t0,tmax,y0,epsilon,h,i):
    t=t0
    n=int((tmax-t0)/h)
    s1=(6,5)
    k=zeros(s1)
    s2=(6,n+1)
    y=zeros(s2)
    yp=zeros(6)
    y[:,0]=y0[:]
    k[:,0]=f(t,y[:,i])
    k[:,1]=f(t+h/3,y[:,i]+h/3*k[:,0])
    k[:,2]=f(t+h/3,y[:,i]+h/6*k[:,0]+h/6*k[:,1])
    k[:,3]=f(t+h/2,y[:,i]+h/2*k[:,0]+h/8*3*k[:,1])
    k[:,4]=f(t+h/2,y[:,i]+h/2*k[:,0]-3/2*h*k[:,2]-2*h*k[:,3]) 
    yp=y[:,i]+h/2*(k[:,0]+3*k[:,2]+4*k[:,3])
    y[:,i+1]=y[:,i]+h/6*(k[:,0]+4*k[:,3]+k[:,4])
    t=t+h
    R=0.2*math.fabs(y[2,i+1]-yp[2])       
    if R>=epsilon:
        i=i+1
        h=h/2.0
        rkm45(func,t0,tmax,y0,epsilon,h,i)
    else:
        if R<=epsilon/64.0:
            i=i+1            
            h=2*h
            rkm45(func,t0,tmax,y0,epsilon,h,i)
        else: 
            while(t<tmax):
                t=t+h                
                k[:,0]=f(t,y[:,i])
                k[:,1]=f(t+h/3,y[:,i]+h/3*k[:,0])
                k[:,2]=f(t+h/3,y[:,i]+h/6*k[:,0]+h/6*k[:,1])
                k[:,3]=f(t+h/2,y[:,i]+h/2*k[:,0]+h/8*3*k[:,1])
                k[:,4]=f(t+h/2,y[:,i]+h/2*k[:,0]-3/2*h*k[:,2]-2*h*k[:,3]) 
                y[:,i+1]=y[:,i]+h/6*(k[:,0]+4*k[:,3]+k[:,4])
                i=i+1
    print y
    return y

def func(t,y):
    nu=zeros(6)
    nu=y
    dk1=nu[1]
    dk3=nu[3]
    dk5=nu[5]
    dk2=F*math.cos(nu[4])/m
    dk4=F*math.sin(nu[4])/m
    dk6=(3*a*F)/((a^2+b^2)*m)
    res=zeros(5)
    res=[dk1,dk2,dk3,dk4,dk5,dk6]
    return res
    


nice=rkm45(func,t0,tmax,y0,epsilon,h,0)



def graf1():
    x = arange(len(nice[0,:]))
    y = nice[0,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()

def graf2():
    x = arange(len(nice[0,:]))
    y = nice[1,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()

def graf3():
    x = arange(len(nice[0,:]))
    y = nice[2,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()

def graf4():
    x = arange(len(nice[0,:]))
    y = nice[3,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()

def graf5():
    x = arange(len(nice[0,:]))
    y = nice[4,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()

def graf6():
    x = arange(len(nice[0,:]))
    y = nice[5,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()
    
graf1()
graf2()
graf3()
graf4()
graf5()
graf6()