蒙特卡洛求定积及龙哥库塔解微分方程
- 格式:doc
- 大小:148.50 KB
- 文档页数:4
龙格库塔法解微分方程组引言微分方程组是数学中经常遇到的问题,在物理、工程和自然科学中都有广泛的应用。
为了求解微分方程组,我们需要利用数值方法来逼近解析解。
本文将介绍一种常用的数值方法——龙格库塔法(Runge-Kutta method),并探讨如何利用该方法来解微分方程组。
龙格库塔法概述龙格库塔法是一种迭代数值方法,用于求解常微分方程的初值问题。
它的主要思想是将微分方程的解进行离散化,将其转化为一系列的逼近值。
龙格库塔法的基本步骤如下:1.确定步长h和迭代次数n。
2.初始化初始条件,并假设第一个逼近值为y(xi)。
3.依次计算每个逼近值,直到得到y(xi+n*h)为止。
在每个迭代步骤中,龙格库塔法根据前一步的逼近值来计算下一步的逼近值。
该方法具有高阶精度和较好的稳定性,在实际应用中广泛使用。
单一微分方程的龙格库塔法首先,我们来看如何利用龙格库塔法来解一阶常微分方程。
以方程dy/dx = f(x, y)为例,其中f(x, y)为给定的函数。
步骤一:确定步长和迭代次数选择合适的步长h和迭代次数n来进行数值计算。
步长h决定了离散化的精度,而迭代次数n决定了逼近解的数目。
步骤二:初始化条件并计算逼近值设初始条件为y(x0) = y0,其中x0为起始点,y0为起始点处的函数值。
我们先通过欧拉法计算出y(x0 + h)的逼近值,然后再通过该逼近值来计算下一个逼近值。
逼近值的计算公式如下:k1 = h * f(x0, y0)k2 = h * f(x0 + h/2, y0 + k1/2)k3 = h * f(x0 + h/2, y0 + k2/2)k4 = h * f(x0 + h, y0 + k3)y(x0 + h) = y0 + 1/6 * (k1 + 2k2 + 2k3 + k4)步骤三:重复步骤二直到得到y(xi+n*h)依次利用上一步计算出的逼近值来计算下一个逼近值,直到得到y(xi+n*h)为止。
微分方程组的龙格库塔法对于一阶微分方程组的初值问题,我们可以将其转化为向量形式。
4阶经典龙格库塔公式求解微分方程4阶龙格库塔法(也称为四阶Runge-Kutta方法)是一种常用于求解常微分方程的数值方法。
它是由德国数学家卡尔·龙格以及他的学生马丁·威尔海姆·库塔于1901年独立提出的。
以下是详细的介绍。
1.问题描述我们考虑一个典型的常微分方程初值问题:dy/dx = f(x, y),y(x0) = y0其中,x0是给定的初始点,y(x)是我们要求解的函数,f(x,y)是已知的函数。
2.方法原理四阶龙格库塔方法的基本思想是通过使用全区间的函数信息来估计下一步的函数值。
在每个步骤中,我们计算四个不同的估计量,然后将它们组合起来以获取最终的解。
具体来说,我们首先计算在初始点x0上的斜率k1=f(x0,y0)。
然后我们计算在x0+h/2处的斜率k2=f(x0+h/2,y0+h/2*k1),其中h是步长。
以此类推,我们计算k3和k4分别在x0+h/2和x0+h处的斜率。
然后,我们通过加权组合这些斜率来计算下一个点的函数值y1:y1=y0+(h/6)*(k1+2*k2+2*k3+k4)。
3.算法步骤以下是使用四阶龙格库塔法求解常微分方程的详细步骤:步骤1:给定初始条件 y(x0) = y0,选择步长 h 和积分终点 x_end。
步骤2:计算积分步数n:n = (x_end - x0)/h。
步骤3:设置变量x=x0和y=y0,并将步长设置为h。
步骤4:对每个步数i=1,2,...,n,执行以下计算:4.1:计算斜率k1=f(x,y)。
4.2:计算斜率k2=f(x+h/2,y+h/2*k1)。
4.3:计算斜率k3=f(x+h/2,y+h/2*k2)。
4.4:计算斜率k4=f(x+h,y+h*k3)。
4.5:计算下一个点的函数值y1=y+(h/6)*(k1+2*k2+2*k3+k4)。
4.6:将x更新为x+h,y更新为y1步骤5:重复步骤4,直到达到积分终点 x_end。
第四讲龙格-库塔⽅法龙格-库塔⽅法3.2 Runge-Kutta法3.2.1 显式Runge-Kutta法的⼀般形式上节已给出与初值问题(1.2.1)等价的积分形式(3.2.1)只要对右端积分⽤不同的数值求积公式近似就可得到不同的求解初值问题(1.2.1)的数值⽅法,若⽤显式单步法(3.2.2)当,即数值求积⽤左矩形公式,它就是Euler法(3.1.2),⽅法只有⼀阶精度,若取(3.2.3)就是改进Euler法,这时数值求积公式是梯形公式的⼀种近似,计算时要⽤⼆个右端函数f的值,但⽅法是⼆阶精度的.若要得到更⾼阶的公式,则求积分时必须⽤更多的f值,根据数值积分公式,可将(3.2.1)右端积分表⽰为注意,右端f中还不能直接得到,需要像改进Euler法(3.1.11)⼀样,⽤前⾯已算得的f值表⽰为(3.2.3),⼀般情况可将(3.2.2)的表⽰为(3.2.4)其中这⾥均为待定常数,公式(3.2.2),(3.2.4)称为r级的显式Runge-Kutta法,简称R-K⽅法.它每步计算r个f值(即),⽽k由前⾯(i-1)个已算出的表⽰,故公式是显式的.例i如当r=2时,公式可表⽰为(3.2.5) 其中.改进Euler 法(3.1.11)就是⼀个⼆级显式R-K ⽅法.参数取不同的值,可得到不同公式.3.2.2 ⼆、三级显式R-K ⽅法对r=2的显式R-K ⽅法(3.2.5),要求选择参数,使公式的精度阶p 尽量⾼,由局部截断误差定义11122211()()[(,())(,)]n n n n n n n T y x y x h c f x y x c f x a h y b hk ++=--+++ (3.2.6)令,对(3.2.6)式在处按Taylor 公式展开,由于将上述结果代⼊(3.2.6)得要使公式(3.2.5)具有的阶p=2,即,必须(3.2.7)即由此三式求的解不唯⼀.因r=2,由(3.2.5)式可知,于是有解(3.2.8)它表明使(3.2.5)具有⼆阶的⽅法很多,只要都可得到⼆阶精度R-K⽅法.若取,则,则得改进Euler法(3.1.11),若取,则得,此时(3.2.5)为(3.2.9)其中称为中点公式.改进的Euler法(3.1.11)及中点公式(3.2.9)是两个常⽤的⼆级R-K⽅法,注意⼆级R-K⽅法只能达到⼆阶,⽽不可能达到三阶.因为r=2只有4个参数,要达到p=3则在(3.2.6)的展开式中要增加3项,即增加三个⽅程,加上(3.2.7)的三个⽅程,共计六个⽅程求4个待定参数,验证得出是⽆解的.当然r=2,p=2的R-K⽅法(3.2.5)当取其他数时,也可得到其他公式,但系数较复杂,⼀般不再给出.对r=3的情形,要计算三个k值,即其中将按⼆元函数在处按Taylor公式展开,然后代⼊局部截断误差表达式,可得可得三阶⽅法,其系数共有8个,所应满⾜的⽅程为(3.2.10)这是8个未知数6个⽅程的⽅程组,解也是不唯⼀的,通常.⼀种常见的三级三阶R-K⽅法是下⾯的三级Kutta⽅法:(3.2.11)附:R-K 的三级Kutta ⽅法程序如下function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满⾜c1+c2+c3=1,(1/6 4/6 1/6)endformat short; 3.2.3 四阶R-K ⽅法及步长的⾃动选择利⽤⼆元函数Taylor 展开式可以确定(3.2.4)中r=4,p=4的R-K ⽅法,其迭代公式为111223344()n n y y h c k c k c k c k +=++++其中1(,)n n k f x y =,2221(,(,))n n n n k f x a h y b hf x y =++,⽽33311322(,)n n k f x a h y b hk b hk =+++ 44411422433(,)n n k f x a h y b hk b hk b hk =++++共计13个参数待定,Taylor 展开分析局部截断误差,使得精度达到四阶,即误差为5()O h 。
龙格库塔法解常微分方程龙格库塔法是一种常用的数值解决常微分方程(ODE)的数值计算方法。
它于1960年末由古斯塔夫·龙格库塔(Gustav Runge-Kutta)提出,并多次得到不同的改进和拓展,成为经典的数值解决ODE的一种方法。
龙格库塔法的基本思想是,将ODE的微分方程化为一组非线性的简单方程,通过求解这些简单方程来近似解决ODE。
龙格库塔法比较通用,可以解决一阶和多阶的常微分方程,其代表性的多阶龙格库塔方程有隐式四阶龙格库塔方程和显式四阶龙格库塔方程。
设微分方程组如下:$$\frac{ d y}{ d x }=f(x,y(x))$$由此可以推导出隐式四阶龙格库塔方程(前向差分形式):$$y_{n+1}=y_{n}+\frac{h}{6}(k_{1}+2k_{2}+2k_{3}+k_{4}) $$其中$h$是步长,$k_{1}=f(x_{n},y_{n})$,$k_{2}=f(x_{n}+\frac{h}{2},y_{n}+\frac{h}{2}k_{1})$,$k_{3}=f(x_{n}+\frac{h}{2},y_{n}+\frac{h}{2}k_{2})$,$k_{4}=f(x_{n}+h,y_{n}+hk_{3})$。
龙格库塔法有一个重要特点,就是具有“四步”方案,即在一个步长下只计算四次,完成整个迭代计算过程。
这使得龙格库塔法可以获得高效、稳定可靠的数值解。
补充来说,不同的函数问题有不同的龙格库塔方程,比如一阶非线性龙格库塔方程、二阶线性龙格库塔方程等。
总的说来,龙格库塔法是一种用来数值解常微分方程的有效、可靠的方法,它可以有效解决一阶、多阶和非线性的微分方程。
以分析性形式构建的方程组都可以得到精确的数值解,而且计算量也比较少,速度较快,是一种常用的数值求解ODE的方法。
因此,龙格库塔法作为一种能够有效解决常微分方程的数值计算方法,在理论和实践中均受到广泛的应用。
相比传统的求解方法,它具有更高的计算效率,能够更快速、更准确地给出ODE的数值解。
常微分方程龙格库塔法在数学的世界里,有一种神秘的生物叫常微分方程。
它们就像是一道道难解的难题,让很多人抓耳挠腮,心里直叫苦。
不过,别担心,今天我们要聊的就是一种解这些难题的法宝——龙格库塔法。
听起来高大上,但其实它并没有那么可怕,反而可以说是我们的好帮手。
想象一下,你在山顶上,俯瞰着山谷。
你能看到小溪、绿树,还有那些蜿蜒的小路。
常微分方程就像是这些小路,虽然看起来复杂,但其实我们只需要找到合适的路径,顺着它一路走下去。
龙格库塔法就像是一双好鞋,能让你在这条路上走得稳稳当当,不用担心摔跤。
你可能会问,什么是龙格库塔法呢?简单来说,它就是一种数值解法,帮助我们在找不到解析解的时候,用一些聪明的方法来近似解决。
这玩意儿有几个版本,最常用的就是四阶龙格库塔法。
你可以把它想象成一个厨师,做菜的时候得先准备好材料,对吧?龙格库塔法也是如此,得先准备好初始条件和方程。
然后,它就开始了它的“烹饪”过程。
先把这些材料混合,取一些小样本,然后再慢慢调味,最后出炉的就是你想要的结果。
想想看,这个过程就像是我们做饭时不断尝味道,直到找到最佳口感。
你可能会觉得,这个方法听起来简单,但它却隐藏着许多智慧。
在每一步中,我们都得计算出一些斜率,这些斜率就像是那条小溪的流速,告诉我们水的流动方向。
通过这些信息,我们就能预测下一个位置在哪里。
每一步都在“拼图”,一点一点把整个图案拼凑起来。
这也挺像我们的生活,逐步向前,调整方向,不断摸索,最终才能看到那幅完整的画面。
这个过程并不是一帆风顺的。
方程可能会“发脾气”,变得特别复杂,让你心里直犯嘀咕。
不过别灰心,龙格库塔法就像是个灵活的解题高手,总能找到突破口。
关键在于,咱们要有耐心,细致入微,才能真正领悟它的奥秘。
数学就像一场旅行,虽然有时会迷路,但只要不放弃,最后总能找到回家的路。
别忘了,随着计算机技术的发展,龙格库塔法也有了更便捷的实现方式。
你只需要轻轻一按,电脑就能帮你完成复杂的计算,简直像是给了你一双“魔法手”。
四阶龙格库塔法求解微分方程组微分方程是数学中的一个重要分支,它描述了自然界中许多现象的发展规律。
在实际应用中,我们经常需要求解微分方程组,以得到系统的演化规律和性质。
本文将介绍一种常用的求解微分方程组的数值方法——四阶龙格库塔法。
一、微分方程组的数值解法微分方程组是形如下式的方程集合:begin{cases} frac{dy_1}{dx}=f_1(x,y_1,y_2,cdots,y_n) frac{dy_2}{dx}=f_2(x,y_1,y_2,cdots,y_n) cdotsfrac{dy_n}{dx}=f_n(x,y_1,y_2,cdots,y_n) end{cases} 其中,$y_1,y_2,cdots,y_n$是关于自变量$x$的未知函数,$f_1,f_2,cdots,f_n$是已知的函数。
求解微分方程组的数值方法主要有以下两种:1.欧拉法欧拉法是最简单的数值方法之一,它的基本思想是将微分方程组中的导数用差分代替,从而得到一个递推公式。
具体而言,欧拉法的递推公式为:$$y_{i+1}=y_i+hf(x_i,y_i)$$其中,$h$是步长,$x_i$和$y_i$分别是第$i$个点的自变量和因变量的值。
欧拉法的精度较低,容易出现数值误差,但是它的计算速度快,适用于一些简单的微分方程组。
2.龙格库塔法龙格库塔法是一种常用的高精度数值方法,它的基本思想是将微分方程组中的导数用一系列加权的差分代替,从而得到一个递推公式。
其中,四阶龙格库塔法是最为常用的一种。
具体而言,四阶龙格库塔法的递推公式为:begin{aligned} k_1&=hf(x_i,y_i)k_2&=hf(x_i+frac{h}{2},y_i+frac{k_1}{2})k_3&=hf(x_i+frac{h}{2},y_i+frac{k_2}{2})k_4&=hf(x_i+h,y_i+k_3)y_{i+1}&=y_i+frac{k_1+2k_2+2k_3+k_4}{6} end{aligned} 其中,$k_1,k_2,k_3,k_4$是加权的差分,$h$是步长,$x_i$和$y_i$分别是第$i$个点的自变量和因变量的值。
龙格-库塔方法8.2.1 二阶龙格-库塔方法常微分方程初值问题:做在点的泰勒展开:这里。
取,就有(8.11)截断可得到近似值的计算公式,即欧拉公式:若取,式(8.11)可写成:或(8.12)截断可得到近似值的计算公式:或上式为二阶方法,一般优于一阶的欧拉公式(8.2),但是在计算时,需要计算在点的值,因此,此法不可取。
龙格-库塔设想用在点和值的线性组合逼近式(8.12)的主体,即用(8.13)逼近得到数值公式:(8.14)或更一般地写成对式(8.13)在点泰勒展开得到:将上式与式(8.12)比较,知当满足时有最好的逼近效果,此时式(8.13)-式(8.14)。
这是4个未知数的3个方程,显然方程组有无数组解。
若取,则有二阶龙格-库塔公式,也称为改进欧拉公式:(8.15)若取,则得另一种形式的二阶龙格-库塔公式,也称中点公式:(8.16)从公式建立过程中可看到,二阶龙格-库塔公式的局部截断误差仍为,是二阶精度的计算公式。
类似地,可建立高阶的龙格-库塔公式,同时可知四阶龙格-库塔公式的局部截断误差为,是四阶精度的计算公式。
欧拉法是低精度的方法,适合于方程的解或其导数有间断的情况以及精度要求不高的情况,当解需要高精度时,必须用高阶的龙格-库塔等方法。
四阶龙格-库塔方法应用面较广,具有自动起步和便于改变步长的优点,但计算量比一般方法略大。
为了保证方法的收敛性,有时需要步长取得较小,因此,不适于解病态方程。
8.2.2 四阶龙格-库塔公式下面列出常用的三阶、四阶龙格-库塔计算公式。
三阶龙格-库塔公式(1)(8.17)(8.18)(3)(8.19)四阶龙格-库塔公式(1)(8.20)(8.21)例8.3用四阶龙格-库塔公式(8.20)解初值问题:解:取步长,计算公式为:计算结果列表8.3中。
表8.3 计算结果1 0.2 1.24789 1.24792 0.000032 3 4 0.40.60.81.637622.296183.533891.637782.296963.538020.000160.000780.004138.2.3 步长的自适应欧拉方法和龙格-库塔方法在计算时仅用到前一步的值,我们称这样的方法为单步法。
华中科技大学硕士学位论文摘要常微分方程广泛出现于物理、生物、医学、工程、控制理论等许多科学与工程领域,其数学表述归结为常微分方程定解问题。
很多偏微分方程问题通过离散空间变量后也可以化为常微分方程问题来近似求解。
因此,研究常微分方程的数值解法具有重要的科学意义。
经过长时间的发展, 常微分方程定解问题的数值解法是已日趋成熟和完善,数值分析工作者构造了许多有实用价值的方法。
本文主要是构造一种新的数值方法,并对其进行理论分析与研究。
在本文的最开始,我们简要介绍Runge-Kutta方法(包括单步Runge-Kutta方法和两步Runge-Kutta方法)的背景,回顾Runge-Kutta方法在解常微分方程初值问题中的阶条件以及稳定性等理论的发展历程。
在第二章,我们首先介绍单步和两步Runge-Kutta方法阶条件。
随后引入了一类三步Runge-Kutta方法,研究了方法的零稳定性。
重点按Albrecht的A--方法推导新构造的三步Runge-Kutta方法的阶条件。
在第三章,我们针对该类三步Runge-Kutta方法,利用第二章推导出的阶条件,构造该类一簇3阶三步Runge-Kutta方法,用计算机搜索寻找到具有较大稳定区间的方法,画出其稳定区域,并与三步显式Adams方法进行比较。
最后一章,我们对新方法进行数值实验,并验证前面的理论分析结果。
关键词:Runge-Kutta方法;阶条件;树根理论;零稳定;绝对稳定区域;A-方法华中科技大学硕士学位论文AbstractOrdinary differential equations (ODEs) are widely used in the fields of physics, biology, medicine and many other engineering fields, its mathematical formulation attribute to ordinary differential equations for solutions. Many partial differential equations can be translated into the approximate solution of ordinary differential equations through discreting the quantity. There, studying the numerical solution of ordinary differential equations is very important. Now the issue of ordinary differential equations for the numerical solution is more and more mature and complete. The workers of numerical analysis constructed many practical value methods. This paper is to construct and research a new numerical method.At the beginning of this paper, we explain the background of the Runge-Kutta method(including single-step and two-step Runge-Kutta method), review the Order conditions and stability theory of Runge-Kutta method in the initial value problems.In the second chapter, we first introduce the order conditions of Runge-Kutta methods, then construct a new type of three-step Runge-Kutta methods and study its zero-stable. We main derive the order conditions of three-step Runge-Kutta method using the A-method of Albrecht.In the third chapter, we aim at that type of three-step Runge-Kutta methods, make use of the order conditions deduced in chapter II , construct three-step Runge-Kutta method of 3-order, searching the method with the calculator to have large regions of absolute stability, draw the region of stability and carry on a comparison with three-step Adams methods.In the last chapter, we carry on the number experiment to the new method, and verify the anterior result.Keywords: Runge-Kutta method; order conditions; algebraic criteria for order;zero-stable; region of absolute stability; A-method独创性声明本人声明所呈交的学位论文是我个人在导师指导下进行的研究工作及取得的研究成果。
毕业论文开题报告信息与计算科学常微分方程初值问题的Runge-Kutta解法一、选题的背景、意义常微分方程是指只有一个自变量的微分方程。
而且常微分方程在很多学科领域内有着重要的作用,如自动控制、各种电子学装置的设计、弹道的计算、飞机和导弹飞行的稳定性的研究、化学反应过程稳定性的研究等等,这些问题都可以化为求微分方程的解,或者化为研究解的性质的问题。
这些问题都包含某个变量关于另一个变量的变化。
大多数这样的问题需要求解一个初值问题,即求解满足给定初值条件的微分方程[]1。
现代科学、技术、工程中的大量数学模型都可以用常微分方程的初值问题来描述,所以研究常微分方程初值问题的一些性质、解法很有必要。
1.1 历史背景常微分方程在微积分概念出现后即已出现。
从莱布尼兹专门研究用变量变换解决一阶微分方程的求解问题的“求通解”时代,到1841年刘维尔的里卡蒂方程不存在一般初等解和柯西的初值问题的提出,常微分方程转向“求定解”时代。
再到19世纪末天体力学的太阳系稳定性问题眼界需要常微分方程解的大范围性态,从而发展到了“求所有解问题”。
最后到了20世纪六七十年代,常微分方程因为计算机的发展进入“求特殊解”阶段,发现了具有新性质的特殊的解和方程,如混沌(解)、奇异吸引子及孤立子等等[]2。
常微分方程初值问题的数值求解是求近似解的一种经典方法。
其中常微分初值问题的数4-和Runge-Kutta 值求解问题的方法又包括单步法和线性多步法[]3。
单步法主要有欧拉法[]8 5-。
Runge-Kutta法是最重要也最便于应用的单步法。
法[]101.2 国内外研究现状常微分方程在科学和工程中常用于建立数学模型,通常它们没有解析解,而需要数值方法来近似求解。
在科学的计算机化进程中,常微分方程也在不断地精进中,更在不同的领域得到了前所未有的发展和应用。
求解常微分方程的初值问题的数值方法有单步法和多步法。
其中单步法中的Runge-Kutta 方法是目前工程上求常微分方程数值解的基本方法。
蒙特卡洛法求定积分:
整体思路,蒙特卡洛求定积分的主要思想就是通过在取值范围内大量随机数的随机选取对函数进行求值,进而除以所取次数并乘以其区间,即为所积值。
Step 1:
在执行程序前,打开matlab,执行以下操作打开File—>New—>Function,并点击Function,弹出Editor窗口,将其中内容修改为
function [ y ] = f( x )
y=cos(x)+2.0;
end
(如图所示)。
Step 2:
在Editor窗口中点击File—>Save As,保存至所建立的文件夹内,保存名称必须为f.m。
Step 3:
回到Matlab程序中,将Current Folder改为与刚刚Function函数定义的保存路径一致的路径。
Step 4:
在Command Window里输入以下源程序。
源程序:
>> n=10^6; %用来表示精度,表示需要执行次数。
>> y=0; %初始化y=0,因为f(x)=cos(x)+2.0,下面出现y=y+f(x),故尔y用来统计计算出的所有f(x)值的和。
>> count=1; %从1开始计数,显然要执行n次。
>> while count<=n %当执行次数少于n次时,继续执行
x=unifrnd(0,4); %在matlab中,unifrnd是在某个区间内均匀选取实数(可为小数或整
数),表示x取0到4的任意数。
y=y+f(x); %求出到目前为止执行出的所有f(x)值的和,并用y表示
count=count+1; %count+1表示已执行次数再加一,将来通过与n的比对,判断是否执行下一次
end %如果执行次数已达n次,那么结束循环
z=y/n*4 %用y除以执行次数n,那么取得平均值,并用它乘以区间长度4,即可得到z 的近似值
z=7.2430
由于matlab中不能使用θ字符,故用z代替,即可得到θ=7.2395。
在执行过程中,发现每一次执行结果都会不一样,这是因为是随机选取,所以不同次数的计算结果会有偏差,但由于执行次数很多,从概率的角度来讲都是与真实值相近似的值。
用四阶龙格库塔法解微分方程:
解:
一、求微分方程的数值解:
方法一:
1.将方程化为1次,即化为:
在执行程序前,打开matlab,执行以下操作打开File—>New—>Function,并点击Function,弹出Editor窗口,将其中内容修改为
将二阶函数化为一阶过程中,定义如下:
function dy=func(~,y)
dy=zeros(2,1); %初始化,2行一列,即()=()
dy(1)=y(2); %将上述方程组化简得来
dy(2)=-4*y(2)-29*y(1); %将上述方程组化简得来
end %定义结束
回到Matlab程序中,将Current Folder改为与刚刚Function函数定义的保存路径一致的路径。
在Command Window里输入以下源程序。
>> [x,y]=ode45('func',[0 12],[0 15])
绘图
>> y=size(y)
y = 229 2% 数组y包含的元素数
>> reshape(y,458,1)%将含有229*2=458格元素的数组y化为458行1列的数组
>> plot(x,y) %绘图
>> title('数值图')%图名
>> hold on%保持当前图形
>> xlabel('x')%加x坐标
>> ylabel('y')%加y坐标
024681012 -10
-5
5
10
15
数值图
x
y
方法二:
2.最常用的R-K公式——标准4阶R-K公式为:
在editor窗口中打开一个新函数,对龙格库塔函数定义:
function [x,y]=lgkt4j(odefun,xspan,y0,h)
%odefun为微分方程的函数描述,xspan为解区间[x0,xn],y0为初始条件,h为迭代步长
x=xspan(1):h:xspan(2); %定义x为从xspan(1)开始到xspan(2)且步长为h
k=1;%初始化k=1
y(:,1)=y0(:); %初始化y(:,1)
for i=1:length(x)-1 %从第一次循环开始,共循环(length(x)-1)次
K1=feval(odefun,x(k),y(:,k));%feval(函数名,x1,x2,…)
K2=feval(odefun,x(k)+h/2,y(:,k)+h/2*K1);
K3=feval(odefun,x(k)+h/2,y(:,k)+h/2*K2);
K4=feval(odefun,x(k)+h,y(:,k)+h*K3); %以上K1,K2,K3,K4均为龙哥库塔定
义中的式子
y(:,k+1)=y(:,k)+h/6*(K1+2*K2+2*K3+K4); %利用龙哥库塔定义通过y(:,k)的
迭代求y(:,k+1)
k=k+1;
end
end
回到Matlab程序中,将Current Folder改为与刚刚Function函数定义的保存路径一致的路径。
在Command Window里输入以下源程序。
>> f=@(x,y)[y(2);-4*y(1)-29*y(1)];%利用刚刚化简的一次方程
>> [x,y]=lgkt4j(f,[0,12],[0,15],1)%利用定义的龙哥库塔,[0,12]为x
区间,[0,0]为()=(),1为
步长。
二、求微分方程的特征解:
关于求微分方程(组)的解析解命令:
dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
关于此公式的注解:在表达微分方程时,用字母D表示求微分,D2、D3等可以表示求高阶微分,任何D后所跟字母为因变量,自变量可以指定或由系统规则选定为缺省。
步骤:在matlab中直接输入源程序:
源程序:
y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') %直接利用解析函数求特解
y =(3*sin(5*x))/exp(2*x) %结果。