龙格库塔法求微分方程2
python龙格库塔法求解二阶方程
Python是一种高级编程语言，广泛应用于科学计算和工程领域。
而在科学计算中，求解二阶常微分方程是一个常见的问题。
龙格库塔法（Runge-Kutta method）是一种常用的数值求解方法，可以用来求解常微分方程的初值问题。
在Python中，我们可以利用其强大的数值计算库来实现龙格库塔法求解二阶常微分方程。
在本文中，我们将介绍如何使用Python中的龙格库塔法来求解二阶常微分方程。
文章将从以下几个方面展开讲解：1. 二阶常微分方程的基本概念2. 龙格库塔法的原理与公式推导3. Python中如何实现龙格库塔法求解二阶常微分方程4. 一个具体的求解例子及代码实现5. 总结与展望一、二阶常微分方程的基本概念二阶常微分方程是指具有形如y''(t) = f(t, y, y')(t)的形式，其中t是自变量，y是未知函数，f是已知函数。
求解二阶常微分方程的目标是找到一个满足方程的函数y(t)。
通常情况下，需要给出该方程的初值条件，即y(0)和y'(0)，以求得方程的具体解。
二、龙格库塔法的原理与公式推导龙格库塔法是一种数值求解常微分方程初值问题的方法，通过迭代计算来逼近方程的解。
它的基本思想是将微分方程的解按照一定的步长进行逼近，然后利用逼近值来计算下一个逼近值，直到达到所需的精度。
龙格库塔法的一般形式为：k1 = h * f(tn, yn)k2 = h * f(tn + 1/2 * h, yn + 1/2 * k1)k3 = h * f(tn + 1/2 * h, yn + 1/2 * k2)k4 = h * f(tn + h, yn + k3)yn+1 = yn + 1/6 * (k1 + 2*k2 + 2*k3 + k4)其中，h是步长，t是自变量，y是未知函数，f是已知函数，k1、k2、k3、k4是中间变量。
三、Python中如何实现龙格库塔法求解二阶常微分方程在Python中，可以利用其强大的数值计算库来实现龙格库塔法求解二阶常微分方程。
龙格库塔求解二阶常微分方程
龙格库塔求解二阶常微分方程一、前言在数学领域中，常微分方程是一种非常重要的数学工具。
它可以用来描述许多自然现象，如物理学、化学、生物学等。
而龙格库塔法是一种求解常微分方程的方法之一。
本文将介绍如何使用龙格库塔法求解二阶常微分方程。
二、什么是二阶常微分方程二阶常微分方程是指形如y''+p(t)y'+q(t)y=f(t)的方程，其中y表示未知函数，p(t)、q(t)和f(t)都是已知函数。
通常情况下，我们需要找到一个特定的y函数来满足这个方程。
三、什么是龙格库塔法龙格库塔法是一种数值求解常微分方程的方法。
它基于欧拉方法，但更准确和稳定。
欧拉方法使用线性逼近来估计未知函数值，而龙格库塔法使用高阶多项式逼近来更准确地估计未知函数值。
四、龙格库塔法的基本思想1. 将时间区间[0,T]平均分成N个小区间；2. 在每个小区间内进行递推计算，直到得到所有时间点上的y值；3. 每次递推计算都使用龙格库塔公式来估算y的值。
五、龙格库塔法的公式对于二阶常微分方程y''+p(t)y'+q(t)y=f(t)，我们可以使用以下公式来递推计算：1. k1=h*y'(t)l1=h*f(t)-p(t)*k1/2-q(t)*y(t);2. k2=h*(y'(t)+l1/2)l2=h*f(t+h/2)-p(t+h/2)*k2/2-q(t+h/2)*(y(t)+k1/2);3. k3=h*(y'(t)+l2/2)l3=h*f(t+h/2)-p(t+h/2)*k3/2-q(t+h/2)*(y(t)+k2/2);4. k4=h*(y'(t)+l3)l4=h*f(t+h)-p(t+h)*k4-q(t+h)*(y(t)+k3);其中，h表示时间步长，t表示当前时间点，f表示已知函数，p和q都是已知常数。
通过这些公式，我们可以得到下一个时间点上的y值。
六、龙格库塔法的实现为了更好地理解龙格库塔法，我们可以看一下具体的实现过程。
二阶微分方程龙格库塔法
二阶微分方程龙格库塔法
1.什么是二阶微分方程？
二阶微分方程是指二阶导数出现的方程，其通用形式为
y''+P(x)y'+Q(x)y=F(x)，其中P(x)、Q(x)、F(x)均为已知函数，y是所求函数。
2.为什么要用龙格库塔法？
求解二阶微分方程是数学、物理等领域中常见的问题，然而大多数情况下无法直接解题，所以需要使用数值方法。
龙格库塔法是一种数值方法，被广泛应用于求解微分方程，其优点是精度高、计算速度快。
3.龙格库塔法的原理
龙格库塔法是一种迭代方法，将微分方程看作初值问题，从初始点出发，采用一定的步长h，在每个点上用插值公式逼近y(x+h)，以此求得y(x+h)的近似值，一步步逼近所要求的精度。
4.龙格库塔法的步骤
（1）确定步长h和积分区间[a,b]。
（2）用初值y(x0)=y0，y'(x0)=y'0求出y(x0+h)的近似值。
（3）根据龙格库塔公式求得y(x0+2h)的近似值。
（4）对于连续求解的情况，重复以上步骤，直到求得所需的精度或者达到指定的终点。
5.龙格库塔公式
龙格库塔法的精度与所采用的公式有关，一般采用二阶或四阶的龙格库塔公式。
二阶龙格库塔公式为：
y0=y(x0)
k1=h*f(x0,y0)
k2=h*f(x0+h,y0+k1)
y(x0+2h)=y0+1/2(k1+k2)
其中，f(x,y)是待求函数。
6.总结
龙格库塔法是求解微分方程的一种常用数值方法，可以高精度、高效地解决二阶微分方程的问题。
该方法所需的计算量较小，易于编写程序实现，在实际应用中具有广泛的用途。
四阶龙格库塔法解微分方程
四阶龙格库塔法解微分方程一、四阶龙格库塔法解一阶微分方程①一阶微分方程：cos y t ，初始值y(0)=0，求解区间[0 10]。
MATLAB 程序：%%%%%%%%%%% 四阶龙哥库塔法解一阶微分方程%%%%%%%%%%% y'=cost%%%%%%%%%%% y(0)=0, 0≤t ≤10，h=0.01%%%%%%%%%%% y=sinth=0.01;hf=10;t=0:h:hf;y=zeros(1,length(t));y(1)=0;F=@(t,y)(cos(t));for i=1:(length(t)-1)k1=F(t(i),y(i));k2=F(t(i)+h/2,y(i)+k1*h/2);k3=F(t(i)+h/2,y(i)+k2*h/2);k4=F(t(i)+h,y(i)+k3*h);y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4)*h;endsubplot(211)plot(t,y,'-')xlabel('t');ylabel('y');title('Approximation');span=[0,10];p=y(1);[t1,y1]=ode45(F,span,p);subplot(212)plot(t1,y1)xlabel('t');ylabel('y');title('Exact');图1②一阶微分方程：()22*/x t x x t =- ，初始值x(1)=2，求解区间[1 3]。
MATLAB 程序： %%%%%%%%%%% 四阶龙哥库塔法解微分方程%%%%%%%%%%% x'(t)=(t*x-x^2)/t^2%%%%%%%%%%% x(1)=2, 1≤t ≤3, h=1/128%%%%%%%%%%% 精确解：x(t)=t/(0.5+lnt)h=1/128; %%%%% 步长tf=3;t=1:h:tf;x=zeros(1,length(t));x(1)=2; %%%%% 初始值F_tx=@(t,x)(t.*x-x.^2)./t.^2;for i=1:(length(t)-1)k_1=F_tx(t(i),x(i));k_2=F_tx(t(i)+0.5*h,x(i)+0.5*h*k_1);k_3=F_tx((t(i)+0.5*h),(x(i)+0.5*h*k_2));k_4=F_tx((t(i)+h),(x(i)+k_3*h));x(i+1)=x(i)+(1/6)*(k_1+2*k_2+2*k_3+k_4)*h; endsubplot(211)plot(t,x,'-');xlabel('t');ylabel('x');legend('Approximation');%%%%%%%%%%%%%%%%%%%%%%%%%%%% ode45求精确解t0=t(1);x0=x(1);xspan=[t0 tf];[x_ode45,y_ode45]=ode45(F_tx,xspan,x0);subplot(212)plot(x_ode45,y_ode45,'--');xlabel('t');ylabel('x');legend('Exact');图2二、四阶龙格库塔法解二阶微分方程①二阶微分方程：cos y t ，初始值y(0)=0，y'(0)=-1，求解区间[0 10]。
四阶龙格库塔解二阶微分方程
四阶龙格库塔解二阶微分方程四阶龙格库塔法是一种常用的数值解微分方程的算法，适用于解一阶、二阶以及高阶微分方程。
本文将围绕“四阶龙格库塔解二阶微分方程”展开，并分步骤进行阐述。
第一步：转换为一阶微分方程组对于二阶微分方程 y''=f(x,y,y')，可以将其转换为一阶微分方程组。
令 u=y'，则 y''=du/dx，原方程变为 du/dx=f(x,y,u)，整理得到以下一阶微分方程组：dy/dx=udu/dx=f(x,y,u)第二步：确定步长和初始条件在使用四阶龙格库塔法时，需要确定一个步长 h，以及初始条件y0 和 u0。
步长 h 决定了每一步积分时横坐标 x 的变化量。
初始条件 y0 和 u0 分别代表了 x=0 时的纵坐标和导数值。
第三步：进行数值积分按照四阶龙格库塔法的步骤进行数值积分，可得到第 k 步的纵坐标 yk 和导数值 uk：k1 = hf(xk,yk,uk)j1 = hukk2 = hf(xk+h/2,yk+j1/2,uk+k1/2)j2 = h(uk+k1/2)k3 = hf(xk+h/2,yk+j2/2,uk+k2/2)j3 = h(uk+k2/2)k4 = hf(xk+h,yk+j3,uk+k3)j4 = h(uk+k3)yk+1 = yk+(j1+2j2+2j3+j4)/6uk+1 = uk+(k1+2k2+2k3+k4)/6其中，k1 到 k4 和 j1 到 j4 分别为由函数 f(x,y,u) 求得的中间变量。
yk 和 uk 即为求得的第 k 步的纵坐标和导数值，yk+1 和uk+1 分别为求得的第 k+1 步的纵坐标和导数值。
第四步：循环迭代将上述计算步骤循环执行至所需的精度或区间内积分完成为止。
即对于每一步 k，根据当前的 yk 和 uk 求解下一步的 yk+1 和 uk+1，并记录下这些变量的值。
至此，关于“四阶龙格库塔解二阶微分方程”的阐述就结束了。
计算方法 15 龙格库塔-常微分方程
以上公式就是三阶龙格-库塔公式
计算方法（2016/2017 第一学期）
西南科技大学
制造科学与工程学院
11
四阶龙格-库塔公式
四阶龙格-库塔公式：
h yn1 yn 6 ( k1 2k2 2k3 k4 ) k1 f ( xn , yn ) h hk1 ) k2 f ( xn , yn 2 2 h hk2 k3 f ( xn 2 , yn 2 ) k4 f ( xn h, yn hk3 )
以上公式就是四阶龙格-库塔公式
也称作经典龙格-库塔公式
计算方法（2016/2017 第一学期） 西南科技大学 制造科学与工程学院
12
高阶龙格-库塔公式
高阶龙格-库塔公式：
yn1 yn h( λ1k1 λ2 k2 λ3 k3 λm km ) k1 f ( xn , yn ) k2 f ( xn 2 h, yn 21hk1 ) k3 f ( xn 3 h, yn 31hk1 32 hk2 ) km f ( xn m h, yn m1hk1 m ( m 1) hkm 1 ) i (i 2,, m) 和 ij (i 2,, m; 其中 i (i 1,, m)，
14
龙格-库塔公式
由于龙格-库塔法的导出基于泰勒展开，故精度
主要受解函数的光滑性影响。
对于光滑性不好的解，最好采用低阶算法来求解，
而将步长 h 取小。
计算方法（2016/2017 第一学期）
西南科技大学
制造科学与工程学院
15
例题
例：利用四阶龙格-库塔方法求解微分方程的
matlab迭龙格库塔法解常微分方程
一、介绍迭龙格-库塔法（Runge-Kutta method）是一种数值求解常微分方程（ODE）的常用方法。
它是由卡尔·迭龙格（Carl Runge）和马丁·威尔黑尔姆·库塔（Wilhelm Kutta）在20世纪初提出的，该方法以两位数值分析家的名字来命名。
二、简单描述迭龙格-库塔法是通过数值逼近的方式，来计算常微分方程的近似解。
它是一种显式求解方法，适用于解非线性常微分方程和具有较大阶数的常微分方程。
三、数学原理迭龙格-库塔法主要是通过将微分方程转化为差分方程，利用数值解的方式来逼近微分方程的解。
它是一种显式方法，通过不断迭代得到下一个时间步的近似解。
四、matlab中的应用在matlab中，可以使用ode45函数来调用迭龙格-库塔法求解常微分方程。
ode45函数是matlab中集成的一个函数，通过调用ode45函数，可以直接求解常微分方程的数值解。
五、实例演示下面通过一个简单的例子来演示如何使用matlab中的ode45函数来求解常微分方程。
我们考虑一个简单的一阶常微分方程：dy/dt = -y初始条件为y(0) = 1。
在matlab中，可以通过以下代码来求解该微分方程：```定义微分方程的函数function dydt = myode(t, y)dydt = -y;调用ode45函数求解[t, y] = ode45(myode, [0, 5], 1);plot(t, y);```运行以上代码，即可得到微分方程的数值解，并通过绘图来展示解的变化。
六、总结迭龙格-库塔法是一种常用的数值解常微分方程的方法，它在matlab中有较为方便的调用方式。
通过ode45函数，可以快速求解常微分方程的数值解，并通过绘图来展示结果。
希望本篇文章对读者有所帮助，谢谢阅读。
七、应用场景和优势在实际应用中，迭龙格-库塔法广泛应用于各种科学和工程领域，如物理学、化学、生物学、经济学等。
四阶龙格库塔方法求解n自由度二阶微分方程
四阶龙格库塔方法求解n自由度二阶微分方程在数值计算中，龙格-库塔方法（Runge-Kutta method）是一种常用的求解常微分方程（ODE）的数值方法。
它是一种多步法（multiple-step method），通过使用不同的近似值来迭代计算ODE 的解。
其中，四阶龙格-库塔方法是最常用的龙格-库塔方法之一，也是一种高阶方法，可用于求解n自由度二阶微分方程。
下面将介绍四阶龙格-库塔方法的原理、步骤和应用。
首先，我们来回顾一下二阶微分方程的一般形式：y''(y)=y(y,y(y),y'(y))这里，y(y,y(y),y'(y))是已知的函数，y(y)是我们要求解的未知函数。
为了求解这个二阶微分方程，我们需要将其转化为一个一阶微分方程的求解问题。
令y(y)=y(y)y(y)=y'(y)这样，我们可以将原始的二阶微分方程转化为如下的一阶微分方程组：y'(y)=y(y)y'(y)=y(y,y(y),y(y))现在，我们可以利用四阶龙格-库塔方法来求解这个一阶微分方程组。
四阶龙格-库塔方法基于泰勒展开的思想，通过使用一系列的近似值来逼近方程的解。
该方法的一般形式为：y1=ℎy(yy,yy)y2=ℎy(yy+ℎ/2,yy+y1/2)y3=ℎy(yy+ℎ/2,yy+y2/2)y4=ℎy(yy+ℎ,yy+y3)yy+1=yy+1/6(y1+2y2+2y3+y4)yy+1=yy+ℎ其中，yy是当前时间步，yy+1是下一个时间步，yy是当前位置的近似解，yy+1是下一个位置的近似解，ℎ是时间步长，y1、y2、y3和y4是中间的近似值。
四阶龙格-库塔方法的步骤如下：Step 1:给定初始条件我们需要提供初始条件，包括初始位置y0和初始速度y0，以及时间步长ℎ、求解的时间区间[y0,yy]和离散时间步数y。
y0=y(y0)y0=y'(y0)Step 2:进行迭代计算从n=0开始，进行y步的迭代计算，其中y=(yy−y0)/ℎ。
(完整版)第二节龙格-库塔方法
因为单步法是 p阶的：h0 , 0 h h0 满足| Tn1 | Chp1
| en1 || en | hL | en | Ch p1 | en |
其中 1 hL, Ch p1
en O(hp )
即取
K * 1K1 2 K2 yi1 yi h(1K1 2 K2 )
其中1 和 2 是待定常数。若取 K1 f ( xi , yi ) ，则
问题在于如何确定 xi p 处的斜率 K2 和常数 1 和 2 。
仿照改进的欧拉方法，用欧拉方法预测 y( xi p ) 的值，
yi p yi phK1
k2
f ( xn
h 2
,
yn
h 2 k1)
hh k3 f ( xn 2 , yn 2 k2 )
k4 f ( xn h, yn hk3 )
例2：用经典的龙格-库塔方法
求解下列初值问题 h 0.1
dy 。 dx
y
2x y
x (0,1)
解：经典的四阶龙格-库塔公式： y(0) 1
E(h) 1 h
绝对稳定域： 1 h 1
当 R 时， 1 h 1 2 h 0
绝对稳定区间：(2, 0)
❖经典的R-K公式：yn1
yn
h 6
(k1
2k2
2k3
k4 )
k1 f ( xn , yn ) yn
k2
f
(
yn1 yn hf ( xn , yn )
n1 yn1 yn1
n1 n h[ f ( xn , yn ) f ( xn , yn )] [1 hf y ( xn , )]n
四阶龙格库塔公式推导过程
四阶龙格库塔公式推导过程好啦，今天咱们来聊聊四阶龙格库塔公式。
这可不是个什么神秘的术语，听上去高大上，其实就是一种数值方法，能帮助咱们解决那些不好解的微分方程。
简单来说，就是给你一把钥匙，帮你打开那些复杂数学问题的大门。
你想想，这就像是把一道难吃的菜变得美味可口，真的特别有成就感。
龙格库塔，这个名字听着就挺酷的，像个超级英雄似的。
不过，它可不是在打怪，而是在帮你数值计算。
四阶龙格库塔，顾名思义，就是它有四个“步骤”。
这些步骤就像是在画一幅画，每一笔都是关键。
咱们平时解微分方程，往往需要找出函数在某个点的值，这个时候，龙格库塔就派上用场了。
咱们得知道这个方法是怎么运作的。
假设你有一个微分方程，咱们需要在某个初始条件下找出函数的值。
这个时候，我们会先从初始值开始，给它取个名字，比如说y0，嘿嘿，简单明了。
咱们设定一个步长，叫h。
这就像是在游泳的时候，每划一下都得算好距离。
然后，龙格库塔的第一步就开始啦，咱们要算出一个中间值k1，这个k1就代表了函数在当前点的变化率。
大家可以想象，这就像是在观察波浪的起伏，能让你了解水的走势。
咱们进入第二步，计算k2。
这一步可有意思了，k2是用k1算出来的，不仅如此，咱们还得把步长h的一半加进去，像是在沙滩上画出新的轨迹。
想象一下，你正准备冲浪，得在合适的时机做出反应，这一步特别重要。
然后第三步k3又来了，类似于k2，不过这回要再加一点变化，咱们总得对波涛汹涌的海浪有点心理准备，不是吗？咱们得算出k4，哇，这一步就像是给这幅画上色。
k4是根据k3的结果再往前走一步，咱们把所有这些k值结合起来，哗啦一声，就能算出下一个点的值。
把这些值加起来，算个平均，这样一来，咱们的结果就出来啦！是不是感觉就像是在做一道色拉，把各种食材混合在一起，最终呈现出一幅美丽的佳肴？听到这里，可能有的小伙伴会想：“这四阶龙格库塔是不是太复杂了？”其实不然，记住这几个步骤，就像背熟了一首歌，慢慢的，咱们就能驾轻就熟。
数值计算：四阶龙格-库塔法for二阶微分方程
数值计算：四阶龙格-库塔法for⼆阶微分⽅程引⾔考虑存在以下⼆阶偏微分⽅程f2⋅¨X(t)+f1⋅˙X(t)+f0⋅X(t)=F(t)如何使⽤四阶龙格-库塔法求解该微分⽅程？⼀阶微分⽅程的解法⾸先回顾下对于⼀阶微分⽅程的解法，现在有以下⼀阶常系数⾮齐次微分⽅程f1⋅˙X(t)+f0⋅X(t)=F(t)⽅程可以写作˙X(t)=F(t)−f0⋅X(t)f1X n+1−X ndt=F(t)−f0⋅X(t)f1X n+1=X n+dt⋅F(t)−f0⋅X nf1=X n+dt⋅f(t,X n)其中dt为时间步长。
四阶龙格-库塔法公式如下：X n+1=X n+dt6k1+2k2+2k3+k4k1=f t,X nk2=f t+dt2,Xn+dt2⋅k1k3=f t+dt2,Xn+dt2⋅k2k4=f t+dt,X n+dt⋅k3利⽤以上公式即可显式推进求解⼆阶微分⽅程的解法现在考虑⼆阶常系数⾮齐次微分⽅程f2⋅¨X(t)+f1⋅˙X(t)+f0⋅X(t)=F(t)⽅程中出现⼆次导数项¨X(t)，难以处理，所以考虑将⼆次导数项¨X(t)转化为⼀次导数项，⾄此，我们引⼊a(t)=X(t)b(t)=˙a(t)=˙X(t)原⽅程可以写成⼀阶偏微分⽅程组{()()()()(){f 2⋅˙b (t )+f 1⋅b (t )+f 0⋅a (t )=F (t )b (t )=˙a(t )仿照以上⼀阶微分⽅程的RK4解法˙a(t )=b (t )˙b (t )=F (t )−(f 1⋅b (t )+f 0⋅a (t ))f 2˙a n +1=a n +dt ⋅b (t )=a n +dt ⋅f (t ,a n ,b n )˙bn +1=b n +dt ⋅F (t )−(f 1⋅b (t )+f 0⋅a (t ))f 2=b n +dt ⋅g (t ,a n ,b n )˙a n +1=a n +dt ⋅f (t ,a n ,b n )˙bn +1=b n +dt ⋅g (t ,a n ,b n )化简⾄此，便可以仿照之前的RK4⽅法进⾏求解，具体公式为a n +1b n +1=a n +dt6⋅(k 1a +2k 2a +2k 3a +k 4a )b n +dt6⋅(k 1b +2k 2b +2k 3b +k 4b )其中，各个系数求解公式为：k 1a k 1b=f (t ,a n ,b n )g (t ,a n ,b n )k 2a k 2b =f (t +dt2,a n +dt2⋅k 1a ,b n +dt2⋅k 1b )g (t +dt2,a n +dt2⋅k 1a ,b n +dt2⋅k 1b )k 3a k 3b=f (t +dt2,a n +dt2⋅k 2a ,b n +dt2⋅k 2b )g (t +dt2,a n +dt2⋅k 2a ,b n +dt2⋅k 2b )k 4a k 4b=f (t +dt ,a n +dt ⋅k 3a ,b n +dt ⋅k 3b )g (t +dt ,a n +dt ⋅k 3a ,b n +dt ⋅k 3b )⾄此已经完成使⽤四阶龙格-库塔法⼆阶微分⽅程求解过程的梳理实例求解，以下偏微分⽅程：\begin{align} 1 \cdot \ddot{X(t)}+0 \cdot \dot{X(t)} +0 \cdot {X(t)} =cos(t)\\ \dot{X(0)}=0\\ {X(0)} =0 \end{align}理论解：\begin{align} {X(t)} =-cos(t)+1 \end{align}#include"RK_ODE.h"#include<iostream>using namespace std;#define M 1{{{{[][]{[][][][][][][][]MatrixXd F2(double t) {MatrixXd res = MatrixXd::Identity(M, M);return res;}MatrixXd F1(double t) {MatrixXd res = MatrixXd::Zero(M, M);return res;}MatrixXd F0(double t) {MatrixXd res = MatrixXd::Zero(M, M);return res;}MatrixXd F(double t) {MatrixXd res = MatrixXd::Ones(M, 1);res(0, 0) = cos(t);return res;}int main() {RK_ODE ode(M, 0.1);ode.X.fill(0.);ode.dX.fill(0.);for (int i = 0; i < 1000; i++) {ode.Solve(F2, F1, F0, F);ode.WriteMotionAndForce("results1.txt", 0); }//ode.SaveSnapshoot("snapshoot.json");cout << ode._Mode << endl;}Processing math: 88%。
【免费下载】数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程
ylabel('角度');
A=[x,y];
%y(find(y==max(y)))
%Num=(find(y==max(y)))
[y,T]=max(y);
fprintf('角度峰值等于%d',y) %角度的峰值也就是 π
fprintf('\n')
fprintf('周期等于%d',T*h)
legend('欧拉');
稳定性很好，RK4 法是四阶方法，每步的误差是 h5 阶，而总积累误差为 h4 阶。所以比欧 拉稳定。
运行第三个程序：在一幅图中显示欧拉法和 RK4 法，随着截断误差的积累，欧拉法产生了 较大的误差 h=0.01
对全部高中资料试卷电气设备，在安装过程中以及安装结束后进行高中资料试卷调整试验；通电检查所有设备高中资料电试力卷保相护互装作置用调与试相技互术关，系电，力根通保据过护生管高产线中工敷资艺设料高技试中术卷资，配料不置试仅技卷可术要以是求解指，决机对吊组电顶在气层进设配行备置继进不电行规保空范护载高与中带资负料荷试下卷高问总中题体资，配料而置试且时卷可，调保需控障要试各在验类最；管大对路限设习度备题内进到来行位确调。保整在机使管组其路高在敷中正设资常过料工程试况中卷下，安与要全过加，度强并工看且作护尽下关可都于能可管地以路缩正高小常中故工资障作料高；试中对卷资于连料继接试电管卷保口破护处坏进理范行高围整中，核资或对料者定试对值卷某，弯些审扁异核度常与固高校定中对盒资图位料纸置试，.卷保编工护写况层复进防杂行腐设自跨备动接与处地装理线置，弯高尤曲中其半资要径料避标试免高卷错等调误，试高要方中求案资技，料术编试交写5、卷底重电保。要气护管设设装线备备置敷4高、调动设中电试作技资气高，术料课中并3中试、件资且包卷管中料拒含试路调试绝线验敷试卷动槽方设技作、案技术，管以术来架及避等系免多统不项启必方动要式方高，案中为；资解对料决整试高套卷中启突语动然文过停电程机气中。课高因件中此中资，管料电壁试力薄卷高、电中接气资口设料不备试严进卷等行保问调护题试装，工置合作调理并试利且技用进术管行，线过要敷关求设运电技行力术高保。中护线资装缆料置敷试做设卷到原技准则术确：指灵在导活分。。线对对盒于于处调差，试动当过保不程护同中装电高置压中高回资中路料资交试料叉卷试时技卷，术调应问试采题技用，术金作是属为指隔调发板试电进人机行员一隔，变开需压处要器理在组；事在同前发一掌生线握内槽图部内 纸故，资障强料时电、，回设需路备要须制进同造行时厂外切家部断出电习具源题高高电中中源资资，料料线试试缆卷卷敷试切设验除完报从毕告而，与采要相用进关高行技中检术资查资料和料试检，卷测并主处且要理了保。解护现装场置设。备高中资料试卷布置情况与有关高中资料试卷电气系统接线等情况，然后根据规范与规程规定，制定设备调试高中资料试卷方案。
二阶龙格库塔法例题
二阶龙格库塔法例题二阶龙格-库塔法（RK2）是一种数值求解常微分方程（ODE）的方法，它基于欧拉法的改进。
下面我将以一个例题来说明二阶龙格-库塔法的应用。
假设我们要求解如下的一阶常微分方程：dy/dx = x + y.初始条件为 y(0) = 1。
我们希望使用二阶龙格-库塔法来计算在 x = 0.1 处的函数值。
首先，我们需要将方程转化为标准形式，即将一阶微分方程表示为 dy/dx = f(x, y) 的形式。
在这个例子中，f(x, y) = x + y。
接下来，我们将计算步长 h，这是我们在每个步骤中前进的距离。
在二阶龙格-库塔法中，我们需要选择一个合适的步长。
通常情况下，我们可以通过试验不同的步长来选择一个合适的值。
在这个例子中，我们选择 h = 0.1。
然后，我们可以开始应用二阶龙格-库塔法的迭代公式。
迭代公式如下：k1 = h f(x, y)。
k2 = h f(x + h/2, y + k1/2)。
y_new = y + k2。
现在，我们可以开始迭代计算。
根据初始条件，我们有 x = 0，y = 1。
我们将进行以下计算步骤：1. 计算 k1 = h f(x, y) = 0.1 (0 + 1) = 0.1。
2. 计算 k2 = h f(x + h/2, y + k1/2) = 0.1 (0.05 + 1 + 0.1/2) = 0.105。
3. 计算 y_new = y + k2 = 1 + 0.105 = 1.105。
这样，我们得到了在 x = 0.1 处的函数值 y_new = 1.105。
通过不断重复上述步骤，我们可以计算出在其他点的函数值。
只需将当前的 x 和 y_new 作为下一步的初始条件，重复上述计算步骤即可。
需要注意的是，二阶龙格-库塔法是一种数值近似方法，其结果并不完全精确。
因此，在实际应用中，我们需要根据问题的要求和精度需求选择合适的步长和迭代次数。
以上就是使用二阶龙格-库塔法求解一阶常微分方程的例子。
常微分方程龙格-库塔方法
得中点公式 （8.2.4）
h h yn1 yn hf xn ，yn f xn，yn ， 2 2
取c=2/3得Heun公式
h 2 2 yn1 yn f xn，yn 3 f xn h，yn hf xn，yn 4 3 3
L
i 1
它与显式R-K公式的区别在于：显式公式中，对系数 aij 求和的上限是 i 1 ，从
aij构成的矩阵是一个严格下三角阵。而在隐式公式中，对系数 aij 求和的上 2， ，L 限是L，从而 aij 构成的矩阵是方阵，需要用迭代法求出近似斜率K i i 1，
而 推导隐式公式的思路和方法与显式公式法类似。
第八章常微分方程数值解法
8.2 Runge-Kutta 方法
8.2.1 Runge-Kutta方法的基本思想
8.2.2 几类显式Runge-Kutta方法
第八章常微分方程数值解法
8.2.1 Runge-Kutta方法的基本思想
显式Euler方法是最简单的单步法，它是一阶的，它可以看作Talylor展开后
f
，而不是高阶导数，将它们作线性
组合得到平均斜率，将其与解的Taylor展开相比较，使前面若干项吻合，从而 得到具有一定阶的方法。这就是 yn1 yn h i K i， Runge-Kutta方法的基本思想，其一般形式 i 1 为 (8.2.1) K1 f xn，yn ，
i 1 Ki f x c h , y c h a K 3， ，L， n i n i ij j ，i 2， j 1 L
到 K i* 。参数i，ci 和 aij 待定，确定它们的原则和方法是：将（8.2.2）式中
* 的 y xn1 在 xn 处作Taylor展开，将 K i 在 xn，y xn 处作二元Taylor展开，
龙格-库塔求解微分方程
MATLAB 在数值分析中的应用摘 要：数值分析是研究数学问题数值解及其理论的一个数学分支,涉及面很广.自计算机在其中应用以后,使其应用价值更为广泛,解决实际问题更为有效.以MA TLAB 为平台的数值分析方法与图形可视化编程在解决一些复杂问题时,大大降低了难度,减少了工作量.本文以经典龙格-库塔算法为核心,介绍用数值方法解常微分方程(组),分析自制系统的稳定性问题.如范德波尔方程.以及混沌现象中奇怪吸引子的模拟实现关键词： MA TLAB; 经典龙格-库塔算法; 解常微分方程(组); 范德波尔方程; 奇怪吸引子;1 MATLAB 平台1.1 矩阵MATLAB 是矩阵实验室（Matrix Laboratory ）的简称, 美国MathWorks 公司出品的商业数学软件.以下介绍在处理方程组时的特殊技法: A=[a b c; d e f; g h i]得到矩阵a b c d e f g h i也可直接输入方程组:function y=equition(t,x)y=[q(2)+q(1)/sqrt(q(1)^2+q(2)^2)*(1-q(1)^2-q(2)^2); % 以q 向量保存方程组;-q(1)+q(2)/sqrt(q(1)^2+q(2)^2)*(1-q(1)^2-q(2)^2)] % q(1)代表应变量x, q(2) 代表应变量y得到1.2 语言MATLAB 编程语言与 C 语言 类似, 格式上略有不同. 里不做详述, 程序中会做必要的标注.但在处理矩阵式十分方便, C++是以数组的形式存放变量,对某个元素都要准确指出其位置.22222222()]()]x y x y x y y x x y x y ⎧=+-+⎪+⎪⎨⎪=-+-+⎪+⎩C++是以数组的形式存放变量,对某个元素都要准确指出其位置,处理循环时算法复杂,即使使用指针也很容易出错.2经典龙格-库塔算法公式(15张)龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。
matlab用四阶龙格库塔法解二阶常微分方程
matlab用四阶龙格库塔法解二阶常微分方程在数学和工程中，常微分方程是描述自然界中各种物理现象和过程的常见数学模型。
常微分方程通常包含未知函数及其导数的方程。
在求解常微分方程时，人们通常使用数值方法来近似求解，其中一种常见的方法是使用龙格-库塔方法。
四阶龙格-库塔方法是一种常用的数值方法，用于求解二阶常微分方程。
它是通过在每个步骤中估计未知函数的导数来数值求解方程。
它的精度比较高，通常被认为是最常用的龙格-库塔方法。
对于一个二阶常微分方程：y''(t) = f(t, y(t), y'(t))其中y(t)是未知函数，f(t, y(t), y'(t))是已知的函数。
我们可以将这个二阶微分方程转化为两个一阶微分方程：y'(t) = u(t)u'(t) = f(t, y(t), u(t))其中u(t)是y'(t)的一个替代函数。
为了使用四阶龙格-库塔方法求解这个方程，我们需要将时间区间[t0, tn]分割为等距的n个子区间，每个子区间的长度为h = (tn - t0)/n。
我们从初始点t0开始，通过多个步骤来逼近解直到tn。
每一步我们使用以下递推公式来估计y(t)和u(t)的值：k1 = h * u(t)l1 = h * f(t, y(t), u(t))k2 = h * (u(t) + l1/2)l2 = h * f(t + h/2, y(t) + k1/2, u(t) + l1/2)k3 = h * (u(t) + l2/2)l3 = h * f(t + h/2, y(t) + k2/2, u(t) + l2/2)k4 = h * (u(t) + l3)l4 = h * f(t + h, y(t) + k3, u(t) + l3)然后我们可以使用以下公式来更新y(t)和u(t)的值：y(t + h) = y(t) + (k1 + 2k2 + 2k3 + k4)/6u(t + h) = u(t) + (l1 + 2l2 + 2l3 + l4)/6这个过程重复n次，直到达到结束时间tn。
数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程
Matlab 应用使用Euler 和Rungkutta 方法解臂状摆的能量方程背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。
由角动量定理我们知道εJ M =化简得到0sin 22=+θθl g dt d 在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为，这样比较容易解。
实际上这是一个解二阶常微分方程的问题。
θ在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上，使用能量法建立方程W T =h mg ∆=2J 21ω化简得到四阶龙格库塔方法使用四级四阶经典显式Rungkutta公式稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。
所以比欧拉稳定。
运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差h=0.01h=0.0001，仍然是开始较为稳定，逐渐误差变大总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。
通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。
三个程序欧拉法clear;clch=0.00001;a=0;b=25;x=a:h:b;y(1)=0;z(1)=0;for i=1:length(x)-1 % 欧拉y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));endplot(x,y,'r*');xlabel('时间');ylabel('角度');A=[x,y];%y(find(y==max(y)))%Num=(find(y==max(y)))[y,T]=max(y);fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h)%周期legend('欧拉');龙格库塔方法先定义函数rightf_sys1.mfunction w=rightf_sys1(x,y,z)w=7.35499*cos(y);clear;clc;%set(0,'RecursionLimit',500)h=0.01;a=0;b=25;x=a:h:b;RK_y(1)=0; %初值%RK_z(1)=0;初值for i=1:length(x)-1K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i));%K1 and L1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);K4=RK_z(i)+h*L3;% K4L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,RK_y,'b+');xlabel('Variable x');ylabel('Variable y');A=[x,RK_y];[y,T]=max(RK_y);legend('RK4方法');fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')%周期fprintf('周期等于%d',T*h)两个方法在一起对比使用跟上一个相同的函数rightf_sys1.mclear;clc; %清屏h=0.0001;a=0;b=25;x=a:h:b;Euler_y(1)=0;%欧拉的初值Euler_z(1)=0;RK_y(1)=0;%龙格库塔初值RK_z(1)=0;for i=1:length(x)-1%先是欧拉法Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h*7.35499*cos(Euler_y(i));%龙格库塔K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 andL1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);% K2 and L2K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);% K3 and L3K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); K4 and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,Euler_y,'r-',x,RK_y,'b-');[y,T]=max(RK_y);%角度的峰值也就是πfprintf('角度峰值等于%d',y)fprintf('\n')%周期fprintf('周期等于%d',T*h)xlabel('时间');ylabel('角度');legend('欧拉','RK4');。
基于龙格库塔求解二阶耦合偏微分方程
基于龙格库塔求解二阶耦合偏微分方程下载提示：该文档是本店铺精心编制而成的，希望大家下载后，能够帮助大家解决实际问题。
文档下载后可定制修改，请根据实际需要进行调整和使用，谢谢!本店铺为大家提供各种类型的实用资料，如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等，想了解不同资料格式和写法，敬请关注!Download tips: This document is carefully compiled by this editor. I hope that after you download it, it can help you solve practical problems. The document can be customized and modified after downloading, please adjust and use it according to actual needs, thank you! In addition, this shop provides you with various types of practical materials, such as educational essays, diary appreciation, sentence excerpts, ancient poems, classic articles, topic composition, work summary, word parsing, copy excerpts, other materials and so on, want to know different data formats and writing methods, please pay attention!1. 引言在科学与工程领域，许多现实问题可以用偏微分方程来建模和描述。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《MATLAB 程序设计实践》课程考核一、编程实现“四阶龙格－库塔（R-K ）方法求常微分方程”，并举一例应用之。
【实例】采用龙格-库塔法求微分方程：⎩⎨⎧==+-=0 , 0)(1'00x x y y y 1、算法说明：在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分方程的初值问题。
其算法公式为：)22(63211k k k hy y n n +++=+其中：⎪⎪⎪⎩⎪⎪⎪⎨⎧++=++=++==) ,()21,21()21 ,21(),(3423121hk y h x f k hk y h x f k hk y h x f k y x f k n n n n n n n n 2、流程图：2.1、四阶龙格－库塔（R-K ）方法流程图：2.2、实例求解流程图：3、源程序代码3.1、四阶龙格－库塔（R-K）方法源程序：function [x,y] = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)%Runge-Kutta 方法解微分方程形为 y'(t)=f(x,y(x))%此程序可解高阶的微分方程。
只要将其形式写为上述微分方程的向量形式%函数 f(x,y): fun%自变量的初值和终值：x0, xt%y0表示函数在x0处的值，输入初值为列向量形式%自变量在[x0,xt]上取的点数：PointNum%varargin为可输入项，可传适当参数给函数f(x,y)%x：所取的点的x值%y：对应点上的函数值if nargin<4 | PointNum<=0PointNum=100;endif nargin<3y0=0;endy(1,:)=y0(:)'; %初值存为行向量形式h=(xt-x0)/(PointNum-1); %计算步长x=x0+[0:(PointNum-1)]'*h; %得x向量值for k=1:(PointNum)%迭代计算f1=h*feval(fun,x(k),y(k,:),varargin{:});f1=f1(:)'; %得公式k1f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:});f2=f2(:)'; %得公式k2f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:});f3=f3(:)'; %得公式k3f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:});f4=f4(:)'; %得公式k4y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6; %得y(n+1)end3.2、实例求解源程序：%运行四阶R-K法clear, clc %清除内存中的变量x0=0;xt=2;Num=100;h=(xt-x0)/(Num-1);x=x0+[0:Num]*h;a=1;yt=1-exp(-a*x); %真值解fun=inline('-y+1','x','y'); %用inline构造函数f(x,y) y0=0; %设定函数初值PointNum=5; %设定取点数[x1,y1]=ode23(fun,[0,2],0);[xr,yr]=MyRunge_Kutta(fun,x0,xt,y0,PointNum);MyRunge_Kutta_x=xr'MyRunge_Kutta_y=yr'plot(x,yt,'k',x1,y1,'b--',xr,yr,'r-')legend('真值','ode23','Rung-Kutta法解',2)hold onplot(x1,y1,'bo',xr,yr,'r*')4、程序运行结果：MyRunge_Kutta_x =0 0.5000 1.0000 1.5000 2.0000MyRunge_Kutta_y =0 0.3932 0.6318 0.7766 0.8645二、编程解决以下科学计算问题：（一）[例7-2-4] 材料力学复杂应力状态的分析——Moore圆。
1、程序说明：利用平面应力状态下斜截面应力的一般公式，画出任意平面应力状态下的应力圆（Moore 圆），求出相应平面应力状态下的主应力（1σ、3σ），并求出该应力状态下任意方位角α的斜截面上的应力σ、τ。
2、程序流程图：σxσyσyσxσατyxτyxτxyτxyτ3、程序代码：clear;clc;Sx=input('Sigma_x(MPa)='); %输入该应力状态下的各应力值Sy=input('Sigma_y(MPa)=');Txy=input('Tau_xy(MPa)=');a=linspace(0,pi,100); %等分半圆周角Sa=(Sx+Sy)/2;Sd=(Sx-Sy)/2;Sigma=Sa+Sd*cos(2*a)-Txy*sin(2*a); %应力圆一般方程Tau=Sd*sin(2*a)+Txy*cos(2*a);plot(Sigma,Tau,Sx,Txy,'.r',Sy,-Txy,'.r'); %画出应力圆，标出该应力状态下各应力参数line([Sx,Sy],[Txy,-Txy]);axis equal; %控制各坐标轴的分度使其相等——使应力圆变圆title('应力圆');xlabel('正应力（MPa）');ylabel('剪应力（MPa）');text(Sx,Txy,'A');text(Sy,-Txy,'B');Smax=max(Sigma);Smin=min(Sigma);Tmax=max(Tau);Tmin=min(Tau);b=axis; %提取坐标轴边界ps_array.Color='k'; %控制坐标轴颜色为黑色line([b(1),b(2)],[0,0],ps_array); %调整坐标轴line([0,0],[b(3),b(4)],ps_array);b=axis; %提取坐标轴边界line([b(1),b(2)],[0,0],ps_array); %画出x坐标轴hold onplot(Sa,0,'.r') %标出圆心text(Sa,0,'O');plot(Smax,0,'.r',Smin,0,'.r',Sa,Tmax,'.r',Sa,Tmin,'.r') %标出最大、最小拉应力、切应力点text(Smax,0,'C');text(Smin,0,'D');text(Sa,Tmax,'E');text(Sa,Tmin,'F');%-----------此部分求某一斜截面上的应力---------------------------------------------t=input('若需求某一截面上的应力，请输入1；若不求应力，请输入0：');while t~=0alpha=input('给出斜截面方向角：alpha=（角度）：');sigma=Sa+Sd*cos(2*(alpha/180*pi))-Txy*sin(2*(alpha/180*pi))tau=Sd*sin(2*(alpha/180*pi))+Txy*cos(2*(alpha/180*pi))plot(sigma,tau,'or')t=input('若还需求其他截面上的应力，请输入1；若要退出，请输入0：');endhold off%----------此部分求该应力状态下的主应力--------------------------------------------Sigma_Max=SmaxSigma_Min=SminTau_Max=TmaxTau_Min=TminSigma1=Smax %得出主应力 Sigma3=Sminl=Sx-Sa;h=Txy;ratio=abs(h/l); %求主应力平面方向角 '主应力平面方向角（角度）：' alpha_0=atan(ratio)/2*180/pi4、程序运行结果：（以MPa 20 MPa, 30 MPa, 100-===xy y x τσσ为例）Sigma_x(MPa)=100 Sigma_y(MPa)=30 Tau_xy(MPa)=-20若需求某一截面上的应力，请输入1；若不求应力，请输入0：1 给出斜截面方向角：alpha=（角度）：30 sigma = 99.8205tau =20.3109若还需求其他截面上的应力，请输入1；若要退出，请输入0：0 Sigma_Max = 105.3087Sigma_Min = 24.6970Tau_Max = 40.3109Tau_Min = -40.2963Sigma1 = 105.3087Sigma3 = 24.6970ans =主应力平面方向角（角度）：alpha_0 = 14.8724（二）实验5（椭圆的交点） 两个椭圆可能具有0 ~ 4个交点，求下列两个椭圆的所有交点坐标：(1) ()()523222=+-+-x y x ; (2) ()()43/3222=+-y x1、算法说明：此题相当于求两个二元二次方程组的解，故为便于清楚地显示出两椭圆的相对位置，用ezplot 函数把两个椭圆画在同一个坐标系中，然后利用fsolve 函数解方程组得到两椭圆的交点即方程组的解。
2、程序流程图：3、程序代码：clear; clc;ezplot('(x-2)^2+(y-3+2*x)^2-5',[-1,5, -8,8]); %画第一个椭圆hold onezplot('2*(x-3)^2+(y/3)^2-4',[-1,5, -8,8]); %画第二个椭圆grid on; %显示网格hold offf1=sym('(x-2)^2+(y-3+2*x)^2=5'); %输入两个椭圆方程f2=sym('2*(x-3)^2+(y/3)^2=4');[x,y]=solve(f1,f2,'x','y'); %联立两个椭圆方程求解交点middle=[x,y]; %合并x,y两个矩阵intersection_x_y=double(middle) %将符号解转换成数值解4、程序运行结果：。