Ռունգե-Կուտտայի մեթոդ

Ռունգե-Կուտտայի մեթոդներ, թվային մեթոդների մեծ դաս` սովորական դիֆերենցիալ հավասարումների և դրանց համակարգերի համար Կոշիի խնդիրը լուծելու համար։ Այս դասի առաջին մեթոդներն առաջարկվել են մոտ 1900 թվականին գերմանացի մաթեմատիկոսներ Կարլ Ռունգեի և Մարտին Վիլհելմ Կուտտայի կողմից։

Ռունգե-Կուտտայի մեթոդների դասին են դասվում Էյլերի ակնհայտ մեթոդը և Էյլերի մոդիֆիկացված մեթոդը վերահաշվարկով, որոնք ներկայացնում են համապատասխանաբար ճշգրտության առաջին և երկրորդ կարգի մեթոդները։ Գոյություն ունեն ճշտության երրորդ կարգի ստանդարտ հստակ մեթոդներ, որոնք չեն ստացել լայն տարածում։ Առավել հաճախ օգտագործվում և տարբեր մաթեմատիկական փաթեթներում (Maple, MathCAD, Maxima) իրագործվել է դասական Ռունգե-Կուտտայի դասական մեթոդը, որն ունի ճշգրտության չորրորդ կարգ։ Բարձր ճշգրտությամբ հաշվարկներ կատարելիս ավելի ու ավելի հաճախ են օգտագործվում ճշգրտության հինգերորդ և վեցերորդ կարգերի մեթոդները։ Ավելի բարձր կարգի սխեմաների կառուցումը լի է հաշվարկային մեծ դժվարություններով։

Յոթերորդ կարգի մեթոդները պետք է ունենան առնվազն ինը փուլ, իսկ ութերորդ կարգի մեթոդները՝ առնվազն 11 փուլ։ 9-րդ և ավելի բարձր կարգերի (չունենալով, սակայն, մեծ գործնական նշանակություն) մեթոդների համար հայտնի չէ, թե որքան փուլեր են անհրաժեշտ ճշտության համապատասխան կարգի հասնելու համար։

Ռունգե-Կուտտայի չորրորդ կարգի դասական մեթոդ 
Ռունգե-Կուտտայի չորրորդ կարգի մեթոդը մշտական ինտեգրման քայլով հաշվարկման ժամանակ այնքան տարածված է, որ այն հաճախ կոչվում է պարզապես Ռունգե-Կուտտայի մեթոդ։

Դիտարկենք Կոշիի խնդիրը սովորական դիֆերենցիալ հավասարումների առաջին կարգի համակարգի համար։ ( հետագայում՝ , а ).

 

Այդ դեպքում հետևյալ կետերում մոտավոր արժեքը հաշվարկվում է իտերացիոն բանաձևով.

 

Նոր արժեքի հաշվարկը տեղի է ունենում չորս փուլով.

 
 
 
 

որտեղ  -ը ցանցի քայլի չափն է -ով։

Այս մեթոդն ունի ճշգրտության չորրորդ կարգ։ Սա նշանակում է, որ սխալը մեկ քայլում ունի   կարգ,  և վերջավոր ինտեգրման միջակայքում ընդհանուր սխալն ունի    կարգը։

Ռունգե-Կուտտայի ակնհայտ մեթոդներ 
Ռունգե-Կուտտայի ակնհայտ մեթոդների ընտանիքը հանդիսանում է ինչպես Էյլերի ակնհայտ մեթոդի, այնպես էլ Ռունգե-Կուտտայի չորրորդ կարգի դասական մեթոդի ընդհանրացումը։ Այն տրվում է հետևյալ բանաձևերով․

 

Որտեղ  h-ն ցանցի քայլի արժեքն է x-ով և նոր արժեքի հաշվարկը տեղի է ունենում s փուլով։

 

Կոնկրետ մեթոդը որոշվում է  թվով  և   և   գործակիցներով։ Այս գործակիցները հաճախ դասակարգվում են աղյուսակում (որը կոչվում է Բուտչերի աղյուսակ)։

 

Ռունգ-Կուտտայի մեթոդի գործակիցները պետք է բավարարեն   պայմանին  համար։ Եթե ​​պահանջվում է, որ մեթոդն ունենա   կարգ, ապա պետք է բավարարի հետևյալ պայմանին․

 

որտեղ  -ը Ռունգե-Կուտտայի մեթոդով ստացված մոտեցումն է։ Բազմակի տարբերակումից հետո այս պայմանը վերածվում է բազմաբնույթ հավասարումների համակարգի՝ մեթոդի գործակիցների վերաբերյալ։

Ռունգե-Կուտտայի ոչ ակնհայտ մեթոդներ 
Մինչ այժմ նշված Ռունգե-Կուտտայի բոլոր մեթոդներն ակնհայտ մեթոդներ են։ Սակայն Ռունգե-Կուտտայի ակնհայտ մեթոդները, որպես կանոն, պիտանի չեն խիստ հավասարումների լուծման համար` նրանց բացարձակ կայունության փոքր շրջանի պատճառով։ Պարզ է, որ Ռունգե-Կուտտայի ակնհայտ մեթոդների անկայունությունը նաև լուրջ խնդիրներ է առաջ բերում մասնակի դիֆերենցիալ հավասարումների ածանցյալների լուծման մեջ։

Ռունգե-Կուտտայի մեթոդների անկայունությունը խթանել է ոչ ակնհայտ մեթոդների զարգացումը։ Ռունգե-Կուտտայի ոչ ակնհայտ մեթոդն ունի հետևյալ տեսքը։

 

որտեղ

 

Հստակ մեթոդը բնութագրական է նրանով, որ   գործակիցների մատրիցան նրա համար ունի ցածր եռանկյունաձև տեսք (ներառյալ զրոյական գլխավոր անկյունագիծը), ի տարբերություն ակնհայտ մեթոդի, որտեղ մատրիցան ունի կամայական ձև։ Դա երևում է նաև Բատչերի աղյուսակում․

 

Այս տարբերության հետևանքն այն է, որ անհրաժեշտ է ամեն քայլի լուծել հավասարումների համակարգը  համար, որտեղ  -ը փուլերի քանակն է։ Սա մեծացնում է հաշվարկային արժեքը, բայց բավարար չափով   դեպքում կարելի է կիրառել սեղմման արտացոլման սկզբունքը և լուծել այս համակարգը պարզ իտերացիայի մեթոդով։ Մեկ իտերացիայի դեպքում սա մեծացնում է հաշվարկային արժեքն ընդամենը երկու անգամ։

Մյուս կողմից, Ժան Կունցմանը (1961) Ջոն Բատչերը (1964) ցույց են տվել, որ ցանկացած  փուլերի ցանկացած քանակի համար գոյություն ունի Ռունգե-Կուտտայի ոչ ակնհայտ մեթոդ ճշգրտության   կարգով։ Սա նշանակում է, օրինակ, որ վերը նկարագրված ակնհայտ քառափուլ չորրորդ կարգի մեթոդի համար գոյություն ունի ոչ ակնհայտ անալոգ ճշգրտության երկու անգամ ավելի մեծ կարգով։

Ռունգե-Կուտտայի երկրորդ կարգի անուղղակի մեթոդ 
Ռունգե-Կուտտայի ամենապարզ ոչ ակնհայտ մեթոդն Էյլերի «վերահաշվարկով» մոդիֆիկացված մեթոդն է։ Այն տրվում է հետևյալ բանաձևով.

Նրա իրականացման համար յուրաքանչյուր փուլում անհրաժեշտ է առնվազն երկու իտերացիա (և ֆունկցիայի երկու հաշվարկ)։

Կանխատեսում.

 .

Ուղղում.

 .

Երկրորդ բանաձևը հավասարեցման համակարգի լուծման պարզ իտերացիա է   համեմատութույամբ, որը գրվել է սեղմված արտացոլման  տեսքով։  Ճշգրտությունը բարձրացնելու համար իտերացիա-ուղղումը կարելի է կատարել մի քանի անգամ՝ փոխարինելով ։  Էյլերի «վերահաշվարկով» մոդիֆիկացված մեթոդն ունի ճշգրտության երկրորդ կարգ։

Հաստատունություն 
Ռունգե-Կուտտայի ոչ ակնհայտ մեթոդների առավելությունն ակնհայտների համեմատությամբ նրանց մեծ հաստատունությունն է, ինչը հատկապես կարևոր է խիստ հավասարումների լուծման ժամանակ։ Որպես օրինակ դիտարկենք y'''  = λy գծային հավասարումը։ Ռունգե-Կուտտայի սովորական մեթոդը, որն օգտագործվում է այս հավասարման համար, հանգում է  , с r իտերացիային, որը հավասար է

 

որտեղ -ն նշանակում միավորների վեկտոր-սյունակ։  ֆունկցիան կոչվում է հաստատունության ֆունկցիա։ Բանաձևը ցույց է տալիս, որ     աստիճանի երկու պոլինոմների հարաբերակցություն է, եթե մեթոդը ունի  փուլ։ Ակնհայտ մեթոդներն ունեն  ցածր եռանկյունաձև խիստ մատրիցա, որտեղից հետևում է, որ , և որ հաստատունության ֆունկցիան բազմանդամ է։

Այս օրինակի թվային լուծումը տալիս է  զրո հետևյալ պայմանի դեպքում՝ , ։ Այդպիսի r–երի բազմությունը կոչվում է բացարձակ հաստատունության ոլորտ։  Մասնավորապես, մեթոդը կոչվում է  A-հաստատուն, եթե բոլոր -երը  գտնվում են բացարձակ հաստատունության ոլորտում։ Ռունգե-Կուտտայի ակնհայտ մեթոդի հաստատունության գործառույթը բազմանդամ է, ուստի Ռունգե-Կուտտայի ակնհայտ մեթոդներն սկզբունքորեն չեն կարող լինել A-հաստատուն։

Եթե ​​մեթոդը  ունի   կարգ, ապա հաստատունության ֆունկցիան բավարարում է   при  պայմանը։  Այսպիսով, հետաքրքրություն է ներկայացնում այս աստիճանի բազմանդամների հարաբերությունը, որը լավագույնս է մոտեցնում էքսպոնենցիալ ֆունկցիան։ Այս հարաբերությունները հայտնի են որպես Պադեի ապրոկսիմացիաներ։  աստիճանի համարիչով  և    աստիճանի  հայտարարով Պադեի ապրոկսիմացիան A հաստատուն է միայն այն ժամանակ, երբ ։

 փուլ ունեցող Հաուս-Լեժանդրի մեթոդն ունի  կարգ,  այդ պատճառով նրա հաստատունության ֆունկցիան Պադեի մոտեցումն է ։  Հետևում է, որ մեթոդը A-հաստատուն է։ Սա ցույց է տալիս, որ   Ռունգե-Կուտտայի A-հաստատուն մեթոդները կարող են լինել կամայականորեն բարձր կարգի։ Ի տարբերություն դրա, Ադամսի մեթոդների A-հաստատունության կարգը չի կարող գերազանցել երկուսը։

 Ծրագրավորման ալգորիթմական լեզուներով լուծման օրինակ 

կատարելով y '= z փոխարինումը և 4y-ը տեղափոխելով աջ կողմ, ստանում ենք հետևյալ համակարգը.

public class MainClass {

	public static void main(String[] args) {
		int k = 2;
		double Xo, Yo, Y1, Zo, Z1;
		double k1, k2, k4, k3, h;
		double q1, q2, q4, q3;
                /*
                 *Նախնական պայմաններ
                 */
		Xo = 0;
		Yo = 0.8;
                Zo = 2;

		h = 0.1; // քայլ

		System.out.println("\tX\t\tY\t\tZ");
		for(; r(Xo,2)<1.0; Xo += h){
			
			k1 = h * f(Xo, Yo, Zo);
			q1 = h * g(Xo, Yo, Zo);
			
			k2 = h * f(Xo + h/2.0, Yo + q1/2.0, Zo + k1/2.0);
			q2 = h * g(Xo + h/2.0, Yo + q1/2.0, Zo + k1/2.0);
			
			k3 = h * f(Xo + h/2.0, Yo + q2/2.0, Zo + k2/2.0);
			q3 = h * g(Xo + h/2.0, Yo + q2/2.0, Zo + k2/2.0);
			
			k4 = h * f(Xo + h, Yo + q3, Zo + k3);
			q4 = h * g(Xo + h, Yo + q3, Zo + k3);
			
			Z1 = Zo + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
			Y1 = Yo + (q1 + 2.0*q2 + 2.0*q3 + q4)/6.0;
			System.out.println("\t" + r(Xo + h, k) + "\t\t" + r(Y1,k) + "\t\t" + r(Z1,k));
			Yo = Y1;
			Zo = Z1;
		}
		
	}
        /**
         * ֆունկցիա կլորացնելու և «պոչը» հեռացնելու համար
         */
	public static double r(double value, int k){
		return (double)Math.round((Math.pow(10, k)*value))/Math.pow(10, k);
	}
        /**
         * ֆունկցիաներ, որոնք ստացվում են համակարգից
         */
	public static double f(double x, double y, double z){
		return (Math.cos(3*x) - 4*y);
	}
	public static double g(double x, double y, double z){
		return (z);
	}

}

using System;
using System.Collections.Generic;

namespace PRJ_RungeKutta
{
    /// <summary>
    ///Ռունգե-Կուտտայի մեթոդի իրականացում սովորական դիֆերենցիալ հավասարման համար
    /// </summary>
    public abstract class RungeKutta
    {
        /// <summary>
        ///Ընթացիկ ժամանակ
        /// </summary>
        public double t;  
        /// <summary>
        /// Պահանջվող որոշում Y [0] - ինքնուրույն որոշում, Y [i] — i-ի ածանցյալ լուծումներ
        /// </summary>
        public double[] Y; 
        /// <summary>
        ///Ներքին փոփոխականներ
        /// </summary>
        double[] YY, Y1, Y2, Y3, Y4;
        protected double[] FY;
        /// <summary>
        /// Конструктор
        /// </summary>
        /// <param name="N">համակարգի չափականություն</param>
        public RungeKutta(uint N)  
        {
            Init(N);
        }
        /// <summary>
        /// Конструктор
        /// </summary>
        public RungeKutta(){}
        /// <summary>
        /// Հիշողության առանձնացում աշխատանքային մասիվների համար
        /// </summary>
        /// <param name="N">Մասիվների չափականություն</param>
        protected void Init(uint N)
        {             
            Y = new double[N];
            YY = new double[N];
            Y1 = new double[N];
            Y2 = new double[N];
            Y3 = new double[N];
            Y4 = new double[N];
            FY = new double[N];
        }
        /// <summary>
        /// Սկզբնական պայմանների սահմանում
        /// </summary>
        /// <param name="t0">Սկզբնական ժամանակ</param>
        /// <param name="Y0">Սկզբնական պայման</param>
        public void SetInit(double t0, double[] Y0) 
        {                                     
            t = t0;
            if (Y == null) 
                Init((uint)Y0.Length);
            for (int i = 0; i < Y.Length; i++)
                Y[i] = Y0[i];
        }
        /// <summary>
        /// Համակարգի աջ մասերի հաշվարկ
        /// </summary>
        /// <param name="t">ընթացիկ ժամանակ</param>
        /// <param name="Y">լուծման վեկտոր</param>
        /// <returns>աջ մաս</returns>
        abstract public double[] F(double t, double[] Y); 
        /// <summary>
        /// Ռունգե-Կուտտա մեթոդի հաջորդ քայլը
        /// </summary>
        /// <param name="dt">ընթացիկ ժամանակի քայլ (կարող է փոփոխական լինել)</param>
        public void NextStep(double dt)
        {
            int i;

            if (dt < 0) return;

            // հաշվարկել Y1
            Y1 = F(t, Y);

            for (i = 0; i < Y.Length; i++)
                YY[i] = Y[i] + Y1[i] * (dt / 2.0);
            
            // հաշվարկել Y2
            Y2 = F(t + dt / 2.0, YY);

            for (i = 0; i < Y.Length; i++)
                YY[i] = Y[i] + Y2[i] * (dt / 2.0);
            
            // հաշվարկել Y3
            Y3 = F(t + dt / 2.0, YY);

            for (i = 0; i < Y.Length; i++)
                YY[i] = Y[i] + Y3[i] * dt;
            
            // հաշվարկել Y4
            Y4 = F(t + dt, YY);

            // հաշվարկել լուծումը նոր քայլով

            for (i = 0; i < Y.Length; i++)
                Y[i] = Y[i] + dt / 6.0 * (Y1[i] + 2.0 * Y2[i] + 2.0 * Y3[i] + Y4[i]);

            //հաշվարկել ընթացիկ ժամանակը
            t = t + dt; 
        }
    }
    class TMyRK : RungeKutta
    {
        public TMyRK(uint N) : base(N) { }

        /// <summary>
        /// օրինակ մաթեմատիկական փարոս 
        /// y''(t)+y(t)=0
        /// </summary>
        /// <param name="t">Ժամանակ</param>
        /// <param name="Y">Լուծում</param>
        /// <returns>Առաջին մաս</returns>
        public override double[] F(double t, double[] Y)
        {
            FY[0] =  Y[1]; 
            FY[1] = -Y[0]; 
            return FY;
        }
        /// <summary>
        ///Օգտագործման օրինակ

        /// </summary>
        static public void Test()
        {
            //Քայլ առ ժամանակ
            double dt = 0.001;
            // Մեթոդի օբյեկտ
            TMyRK task = new TMyRK(2);
            // Սահմանել նախնական պայմանները y (0)=0, y'(0)=1 խնդիր
            double[] Y0 = { 0, 1 }; 
            // Սահմանենք խնդրի սկզբնական պայմանները
            task.SetInit(0, Y0);
            // լուծել մինչև 15 վայրկյան
            while (task.t <= 15) 
            {
                Console.WriteLine("Time = {0:F5}; Func = {1:F8};  d Func / d x = {2:F8}", task.t, task.Y[0], task.Y[1]); // դուրս բերել
t, y, y'
                //հաշվարկել հաջորդ քայլը, քայլ ինտեգրման
                task.NextStep(dt);
            }
            Console.ReadLine();
        }
    }

    class Program
    {
        static void Main(string[] args)
        {
            TMyRK.Test();
        }
    }
}

C # ծրագրում օգտագործվում է RungeKutta աբստրակտ դասը, որի մեջ պետք է վերափոխվի F վերացական մեթոդը, որը սահմանում է հավասարումների աջ կողմերը։

 MATLAB միջավայրում լուծման օրինակ 
Դիֆերենցիալ հավասարումների համակարգերի լուծումը Ռունգե-Կուտտայի մեթոդով հանդիսանում է լուծումների ամենատարածված թվային մեթոդներից մեկը տեխնիկայում։ MATLAB-ի միջավայրում իրականացվել է նրա տեսակներից մեկը՝ Դորման-Փրինսի մեթոդը։ Հավասարումների համակարգի լուծման համար անհրաժեշտ է նախ գրել ածանցյալները հաշվող ֆունկցիաները, այսինքն` y = g(x, y, z) և z = cos(3x) − 4y = f(x, y, z) ֆունկցիաները, ինչի մասին նշված է վերևում։ Կատալոգներից մեկում, որը հասանելի է MATLAB համակարգից, պետք է ստեղծել տեքստային ֆայլ runge.m անունով (օրինակ)  և հետևյալ բովանդակությամբ (MATLAB 5.3 տարբերակի համար)։

Ֆայլի անվանումը և ֆունկցիայի անվանումը պետք է լինեն նույնը, բայց այն կարող է լինել նախկինում չօգտագործված։

Այնուհետև պետք է ստեղծել հիմնական ֆայլ, օրինակ՝ main.m անունով, որը կկատարի հիմնական հաշվարկները։ Այս հիմնական ֆայլը պարունակում է հետևյալ տեքստը՝

Քանի որ MATLAB-ը կենտրոնացած է մատրիցների հետ աշխատելու համար, Ռունգե-Կուտտայի մեթոդի լուծումը շատ հեշտ է կատարվում մի շարք x-ի համար, ինչպես, օրինակ, ծրագրի բերված օրինակում։ Այստեղ լուծումն է՝ ֆունկցիայի գրաֆիկը 0-ից մինչև x_fin ժամանակների սահմանում.

ODE45 ֆունկցիայի արդյունքում ստացված x և y փոփոխականները արժեքային վեկտորներ են։ Ակնհայտ է, որ վերը տրված կոնկրետ օրինակի լուծումը երկրորդ x տարրն է, քանի որ առաջին արժեքը 0 է, ինտեգրման քայլը՝ h = 0.1, իսկ որոնվող արժեքը՝ x = 0.1: MATLAB հրամանի պատուհանում հետևյալ տեքստի մուտքագրումը կտա ցանկալի լուծում.Պատասխան՝ y1 = 0,98768''

Ծանոթագրություններ

Գրականություն