龙格库塔法求微分方程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)/ℎ。
- 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、程序运行结果：。