四阶龙格-库塔(R-K)方法求常微分方程
四阶龙格库塔实验报告
三、四阶Runge-Kutta 法求解常微分方程一、龙格库塔法的思想根据第九章的知识可知道，Euler 方法的局部截断误差是2()O h ，而当用Euler 方法估计出1,()(1)n n n n y y hf x y +=+ 再用梯形公式111[(,)(,)](2)2n n n n n n h y y f x y f x y +++=++进行校正，即采用改进Euler 方法得出数值解的截断误差为3()O h 。
由Lagrange 微分中值定理'11()()()()()(,())(3)n n n n n y x y x y x x y x hf y ξξξ++=+-=+ 记*(,())k hf y ξξ=，得到*1()()(4)n n y x y x k +=+这样只要给出一种计算*k 的算法，就能得到相应的计算公式。
用这种观点的来分析Euler 方法和改进Euler 方法，Euler 方法的迭代公式可改写为111(,)n n n n y y k k hf x y +=+=改进Euler 方法的预报-校正公式可改写为 1121211()2(,),(,)n n n n n n y y k k k hf x y k hf x h y k +=++==++ Euler 方法实际上是用一个点处的值1k 近似*k ，而改进Euler 方法是用两个点处的值1k ，和2k ，做算术平均值近似*k 自然改进Euler 方法要优于Euler 方法。
因此，可以想到假如在1[,]n n x x +内多预报几个点值i k ，并用他们的加权平均值作为*k 的近似值，则有可能构造出具有更高精度的计算公式，这就是Runge-Kutta 法的基本思想。
二、四阶龙格库塔法由Runge-Kutta 的基本思想，构造四阶Runge-Kutta 法是利用1234,,k k k k 和的加权平均值来近似*k ，因此令1112233441211132211243312213(,)(,)(5)(,)(,)n n n n n n n n n n y y w K w K w K w K K hf x y K hf x h y K K hf x h y K K K hf x h y K K K αβαβγαβγη+=++++⎧⎪=⎪⎪=++⎨⎪=+++⎪⎪=++++⎩使得511()()n n y x y O h ++-=即其总体截断误差为4()O h 。
四阶龙格-库塔法求解常微分方程的初值问题-matlab通用程序
参考教材《数值分析》李乃成.梅立泉clearclcformat longm=input('请输入常微分方程的阶数m=');a=input('请输入x下限a=');b=input('请输入x上限b=');h=input('请输入步长h=');ym=input('令y(1,1)=y,y(2,1)=y’,y(3,1)=y’’...请输入ym=','s'); %输入的时候必须按照这个形式输入y1=y(1,1);if m==1 %一阶初值问题单独求解mm=(b-a)/h;y(1,1)=input('请输入在初值点的函数值f(a)=');x=a;y11(1)=y(1,1);for k1=2:(mm+1)y1=y(1,1);K(1,1)=h*(eval(ym)); %计算K1x=x+h/2;y(1,1)=y1+K(1,1)/2;y1=y(1,1);K(1,2)=h*(eval(ym)); %计算K2x=x;y(1,1)=y1+K(1,2)/2-K(1,1)/2;y1=y(1,1);K(1,3)=h*(eval(ym)); %计算K3x=x+h/2;y(1,1)=y1+K(1,3)-K(1,2)/2;y1=y(1,1);K(1,4)=h*(eval(ym)); %计算K4y11(k1)=y11(k1-1)+(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4))/6; y(1,1)=y11(k1);x=a+(k1-1)*h;endy11else %高阶初值问题mm=(b-a)/h; %一共要求解mm个数据点for k2=1:m %读取初值条件fprintf('请输入%d阶导数的初值f(%d)(a)=\n',(k2-1),(k2-1));y(k2,1)=input('=');endfor k2=1:my22(1,k2)=y(k2,1); %先把初值保存在矩阵y22(m,n)中,m表示第几个所求点,n表示第n阶初值endx=a;for k4=2:(mm+1) %求解mm个数据点的循环for k=1:(m-1) %计算K1,包括每一阶的K1 K(k,1)=h*y(k+1,1); %y(k+1,1)中k+1表示第k+1阶,1表示第一个点;K(k,1)中k表示阶数,1表示K1endK(m,1)=h*(eval(ym));x=x+h/2; %求解K1之前,先重新对x和y赋值for k3=1:my(k3,1)=y(k3,1)+K(k3,1)/2;endfor k=1:(m-1) %计算K2K(k,2)=h*y(k+1,1);endK(m,2)=h*(eval(ym));x=x;for k3=1:my(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2;endfor k=1:(m-1) %计算K3K(k,3)=h*y(k+1,1);endK(m,3)=h*(eval(ym));x=x+h/2;for k3=1:my(k3,1)=y(k3,1)+K(k3,3)-K(k3,2)/2; %这里容易出错endfor k=1:(m-1) %计算K4K(k,4)=h*y(k+1,1);endK(m,4)=h*(eval(ym));for k5=1:my22(k4,k5)=y22(k4-1,k5)+(K(k5,1)+2*K(k5,2)+2*K(k5,3)+K(k5,4))/6; %这里,除了要求出下一个点的数值,还要求出相应的导数值endfor k6=1:m %除了对y(1,1)重新赋值外,还要对y(2,1)等重新赋值y(k6,1)=y22(k4,k6);endx=a+(k4-1)*h;endy22(:,1) end。
四阶Runge-Kutta方法
实验题目3 四阶Runge-Kutta 方法摘要一阶常微分方程初值问题0(,)()dy f x y dxy x y ⎧=⎪⎨⎪=⎩ （6.1） 的数值解法是近似计算中很重要的部分。
常微分方程初值问题的数值解法是求方程（6.1）的解在点列1(0,1,)n n n x x h n -=+=L 上的近似值n y ，这里n h 是1n x -到n x 的步长，一般略去下标记为h 。
常微分方程初值问题的数值解法一般分为两大类：（1）单步法：这类方法在计算n y 时，只用到1n x +、n x 和n y ，即前一步的值。
因此，在有了初值以后就可以逐步往下计算。
典型方法如龙格–库塔()R K -方法。
（2）多步法：这类方法在计算1n y +时，除用到1n x +、n x 和n y 以外，还要用(1,2,,;0)n p y p k k -=>L ，即前面k 步的值。
典型方法如Adams 方法。
经典的R K -方法是一个四阶的方法，它的计算公式是：112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ （6.2） R K -方法的优点是：单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值f 。
前言利用四阶龙格-库塔方法求解微分方程的初值问题程序设计流程龙格-库塔法流程图问题1（1）TestRK4('ode1', 1, [0 -1], 5, inline('-x-1'))TestRK4('ode1', 1, [0 -1], 10, inline('-x-1'))TestRK4('ode1', 1, [0 -1], 20, inline('-x-1'))（2）TestRK4('ode2', 1, [0 1], 5, inline('1./(x+1)'))TestRK4('ode2', 1, [0 1], 10, inline('1./(x+1)'))TestRK4('ode2', 1, [0 1], 20, inline('1./(x+1)'))问题2（1）TestRK4('ode3', 3, [1 0], 5, inline('x.^2.*(exp(x)-x)'))TestRK4('ode3', 3, [1 0], 10, inline('x.^2.*(exp(x)-x)'))TestRK4('ode3', 3, [1 0], 20, inline('x.^2.*(exp(x)-x)'))（2）TestRK4('ode4', 3, [1 -2], 5, inline('2*x./(1-2*x)'))TestRK4('ode4', 3, [1 -2], 10, inline('2*x./(1-2*x)'))TestRK4('ode4', 3, [1 -2], 20, inline('2*x./(1-2*x)'))问题3（1）TestRK4('ode5', 1, [0 1/3], 5, inline('x.^2+1/3*exp(-20*x)'))TestRK4('ode5', 1, [0 1/3], 10, inline('x.^2+1/3*exp(-20*x)'))TestRK4('ode5', 1, [0 1/3], 20, inline('x.^2+1/3*exp(-20*x)'))（2）TestRK4('ode6', 1, [0 1], 5, inline('exp(-20*x)+sin(x)'))TestRK4('ode6', 1, [0 1], 10, inline('exp(-20*x)+sin(x)'))TestRK4('ode6', 1, [0 1], 20, inline('exp(-20*x)+sin(x)'))（3）TestRK4('ode7', 1, [0 0], 5, inline('exp(x).*sin(x)'))TestRK4('ode7', 1, [0 0], 10, inline('exp(x).*sin(x)'))TestRK4('ode7', 1, [0 0], 20, inline('exp(x).*sin(x)'))实验所用函数function [x,y] = RK4ODE(fun, xEnd, ini, h)% RK4ODE 用四阶Runge-Kutta法解初值问题dy/dx = f(x,y),y(x0) = y0，在x处y的值%% Synopsis: [x,y] = RK4ODE(fun, xEnd)% [x,y] = RK4ODE(fun, xEnd, ini)% [x,y] = RK4ODE(fun, xEnd, ini, h)%% Input: fun = (string) 初值问题的函数% xEnd = 使用Euler法的截止点% ini = (optional)初始条件[x0 y0]，默认为[0 0]% h = (optional)步长，默认为0.05%% Output: y = 初值问题在x处y的近似值if nargin < 3ini = [0 0]; %若未给初始条件，将初始条件设为[0 0]endif nargin < 4h = 0.05; %若未给出步长，将步长设为0.05endini = ini(:); %将初始条件转为列向量，便于判断是否正确[m,n] = size(ini);if m ~= 2 | n~= 1error('初始值必须是一个含两个元素的向量[x0 y0]');endx0 = ini(1); %初始化x0y0 = ini(2); %初始化y0x = (x0:h:xEnd)'; %构建x向量y = y0*ones(length(x), 1); %初始化y向量for j=2:length(x)k1 = h * feval(fun, x(j-1), y(j-1)); %三阶Runge-Kutta法的递推公式：y(n+1) = y(n) + (k1 + 2*k2 + 2*k3 + k4) / 6k2 = h * feval(fun, x(j-1)+h/2, y(j-1)+k1/2); % k1 = h * f( x(n), y(n) )k3 = h * feval(fun, x(j-1)+h/2, y(j-1)+k2/2); % k2 = h * f( x(n)+h/2, y(n)+k1/2 )k4 = h * feval(fun, x(j-1)+h, y(j-1)+k3); % k3 = h * f( x(n)+h/2, y(n)+k2/2 )y(j) = y(j-1) + (k1+2*k2+2*k3+k4)/6; % k4 = h * f( x(n)+h, y(n)+k3 )endfunction TestRK4(fun, xEnd, ini, N, result)h = (xEnd - ini(1))/N;[x,y] = RK4ODE(fun, xEnd, ini, h);y0 = feval(result, x);plot(x,y,'ro',x,y0,'b-');function dydx = ode1(x,y)dydx = x+y;function dydx = ode2(x,y)dydx = -y.^2;function dydx = ode3(x,y)dydx = 2*y./x + x.^2.*exp(x);function dydx = ode4(x,y)dydx = (y.^2 + y)./x;function dydx = ode5(x,y)dydx = -20*(y-x.^2) + 2*x;function dydx = ode6(x,y)dydx = -20*y + 20*sin(x) + cos(x);function dydx = ode7(x,y)dydx = -20*(y - exp(x).*sin(x)) + exp(x).*(sin(x) + cos(x));思考题1、实验一中(1)数值解和解析解相同，(2)数值解和解析解稍有不同，因为四阶Runge-Kutta方法是以小段的线性算法来近似获得微分方程的数值解，(1)的准确解是1阶的，(2)的准确解是无限阶的，因此对于(1)数值解和解析解相同。
龙格-库塔法,求解常微分方程
隆格库塔法求解常微分方程摘要科学技术中常常需要求解常微分方程的定解问题，这里问题最简单的形式，是本章将着重考察的一阶方程的初值问题.虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法求解.本文着重讨论了隆格库塔法求解一阶常微分方程的初值问题，采用了精度较高的经典的四阶隆格库塔法,然后通过对实例运用Matlab编程进行计算求解，为了体现计算结果的精确性和方法的优越性，再采用了欧拉法和预估较正法对实例进行计算求解作为比较.通过比较三种方法的计算精度，发现四阶经典龙格-库塔方法的误差最小，预估较正法其次，欧拉方法误差则比较大.最后通过选取不同的步长，研究了不同的步长对隆格库塔法求解常微分方程初值问题的计算精度的影响.总之，本文全面分析了隆格库塔法在求解常微分方程的应用，相比与其他的数值解法，隆格库塔法计算精度较高，收敛性较好，其中四阶的隆格库塔法的效率最高，精度也最高.关键词：四阶隆格库塔法；欧拉法；预估较正法；一阶常微分方程；MatlabRunge Kutta Method For Solving Ordinary Differential EquationsABSTRACTProblem solving ordinary differential equations are often needed in science andtechnology. the problem in the simplest form is the initial value problem of first order equations in this chapter ,which will be discussed. Although there are various analytical methods for solving ordinary differential equations, the analytical method can only be used to solve some special types of equations.differential equations can be summed up the actual problems whichThis paper discusses the initial value problem of Runge Kutta Barclays by solving a differential equation, using the four order Runge Kutta method with high accuracy.for instance through classic Matlab programming calculation, the superiority in order to accurately and reflect the calculation result, then the Euler method and the prediction correction method for instance by calculation through the calculation precision. The comparison of three kinds of methods, found that the error of four order Runge Kutta method of minimum, prediction correction method secondly, Euler method error is relatively large. Finally, by selecting different step, study the affect the calculation accuracy of different step of Runge Kutta method to solve initial value problems of ordinary differential equations.In short, this paper comprehensively analyzes the application of Runge Kutta method for solving ordinary differential equations, compared with the numerical solution of other, higher accuracy Runge Kutta method, good convergence, the Runge Kutta method of order four of the highest efficiency and its precision is the highest.Key words: Four order Runge Kutta method; Euler method; prediction correction method; first order ordinary differential equations; Matlab目录1 问题的提出 (1)1.1 问题背景............................................... . (1)1.2 问题的具体内容 (1)2 问题假设 (2)3 符号系统 (2)4 问题的分析 (3)4.1 欧拉格式 (3)4.2 预估较正法 (3)4.3 四阶隆格库塔法的格式 (4)5 模型的建立与求解 (4)5.1 隆格库塔法的基本原理 (4)5.1.1 Taylor级数 (4)5.1.2 隆格库塔法的基本思想 (4)5.1.3 四阶的隆格库塔法 (5)5.2 其他求解常微分方程边值问题算法的简介 (6)5.3 模型求解 (8)5.3.1 运用MATLAB软件对模型求解结果及析 (8)6 模型的评价 (16)7 课程设计的总结与体会 (16)参考文献 (17)附录 (18)一、问题的提出1.1 问题背景：科学技术中常常需要求常微分方程的定解问题，微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：0(,)()dy f x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩ （1）其中a ，b 为常数.虽然求解此类微分方程有各种各样的解析方法，但解析方法只能用于求解一些特殊类型方程，实际问题中归结出来的微分方程主要靠数值解法求解.因为一阶常微分方程简单但又是求解其他方程的基础，所以发展了许多典型的解法.本文着重讨论一类高精度的单步法——隆格库塔法，并且运用四阶的隆格库塔格式来求解初值问题，并且通过实例运用四阶的隆格库塔格式来求解初值问题，同时与显式与隐式的Euler 格式求解出的结果进行精度比较.1.2 问题的具体内容实例一：在区间[0,1]上采用经典的四阶隆格库塔方法求解微分方程1(0)1dy y x dx y ⎧=-++⎪⎨⎪=⎩，其精确解为x y x e -=+，步长为0.5,然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.实例二：在区间[0,1]上用经典的四阶龙格库塔方法求解初值问题 (0)1x dy xe y dx y -⎧=-⎪⎨⎪=⎩， 其精确解为21(2)2xx e -+，然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.最后在区间[0,1]上分别取步长h=0.1;0.05时进行计算，并且探究选取不同的步长对计算结果精度的影响.二、问题假设2.1 假设数值方法本身的计算是准确的.2.2 假设选取的步长趋于0时计算的结果会收敛到微分方程的准确解.2.3 假设步长的增加不会导致舍入误差的严重积累.三、符号系统3.1 符号说明符号含义h选取的步长*K平均斜率p精度的阶数∆前后两次计算结果的偏差y第n个节点的实验值n()y x第n个节点的精确值nδ实验值与精确值的绝对误差四、问题的分析问题要求运用隆格库塔算法来求解一阶微分方程的初值问题，针对前面提出的实例，本文先用经典的四阶隆格库塔法来求解上面的微分方程，为了体现隆格库塔法的优越性，同时用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.最后在区间[0,1]上分别取步长h=0.1;0.05时进行计算，分析在选取不同的步长时，求解结果的精度变化如何.下面是欧拉法，预估校正法以及经典的四阶隆格库塔法的计算公式.4.1欧拉格式（1）显式欧拉格式1(,)n n n n y y hf x y +=+ （2） 局部截断误差22211()()()()22n n n h h y x y y y x o h ξ++''''-=≈= （3） （2）隐式欧拉格式111(,)n n n n y y hf x y +++=+ （4）局部截断误差2211()()()2n n n h y x y y x o h ++''-≈-= (5) 4.2 预估校正法预估 1(,)n n n n y y hf x y +=+ (6)校正 111[(,)(,)]2n n n n n n h y y f x y f x y +++=++ (7) 统一格式 1[(,)(,(,))]2n n n n n n n n h y y f x y f x h y hf x y +=++++ (8) 平均化格式 11(,),(,),1().2p n n n c n n p n p c y y hf x y y y hf x y y y y ++⎧⎪=+⎪=+⎨⎪⎪=+⎩ (9)4.3 四阶龙格库塔方法的格式（经典格式）112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩(10)五、模型的建立与模型求解5.1 隆格库塔法的基本原理隆格库塔法是一种高精度的单步法，这类方法与下述Taylor 级数法有着紧密的联系.5.1.1 Taylor 级数设初值问题 '00(,)()y f x y y x y ⎧=⎨=⎩有解，按泰勒展开，有2'''1()()()()....2n n n n h y x y x hy x y x +=+++； （11） 其中()y x 的各阶导数依据所给方程可以用函数f 来表达，下面引进函数序列(,)j f x y 来描叙求导过程，即(0)(1)'(0)''(1)(1)(1)'''(2),f f y f f y f f x x f f y f f x y ∂∂=≡=+≡∂∂∂∂=+≡∂∂ （12）(2)(2)()(1)j j j j f f y f f x y ---∂∂=+≡∂∂ （13） 根据上式，结果导出下面 Taylor 格式2'''()1...2!!pp n n nn n h h y y hy y y p +=++++ （14）其中一阶Taylor 格式为： '1n n n y y hy +=+提高Taylor 格式的阶数p 即可提高计算结果的精度，显然，p 阶Taylor 格式的局部截断误差为：11(1)1(1)!p p n n h y y y p ζ+++-+=+ （15） 因此它具有p 阶精度.5.1.2 隆格库塔方法的基本思想隆格库塔法实质就是间接地使用Taylor 级数法的一种方法,考察差商1()()n n y x y x h+- 根据微分中值定理，这有'1()()()n n n y x y x y x h h θ+-=+ （16） 利用所给方程 '(,)y f x y =1()()(,())n n n n y x y x hf x h y x h θθ+=+++ （17） 设 平均斜率*(,())n n K f x h y x h θθ=++，由此可见，只要对平均斜率一种算法，便相应地可以导出一种计算格式.再考察改进的Euler 格式，它可以改写成平均化的形式：1121211()2(,)(,)n n n n n n h y y K K K f x y K f x y hK ++⎧=++⎪⎪=⎨⎪=+⎪⎩（18） 这个处理过程启示我们，如果设法在1(,)n n x x +内多预测几个点的斜率值，然后将它们加权平均作为平均斜率，则有可能构造具有更高精度的计算格式，这就是隆格库塔法的基本思想.5.1.3 四阶的隆格库塔法为了方便起见，本文主要运用经典的隆格库塔算法-四阶隆格库塔格式.其格式如下：112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩（19）下面为其具体的算法流程图：5.2 其他求解常微分边值问题算法的简介5.2.1欧拉数值算法（显式）微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：(,)()dyf x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩(20) 其中a ，b 为常数.开始输入x 0,y 0,h,N x 1=x 0+hk 1=f(x 0,y 0),k 2=f(x 0+h/2,y 0+hk 1/2)k 3=f(x 0+h/2,y 0+hk 2/2),k 4=f(x 1,y 0+hk 3)y 1=y 0+h(k 1+2k 2+2k 3+k 4)/6n=1输出x 1,y 1n=N? 结束n=n+1x 0=x 1,y 0=y 1否是图5.1 龙格-库塔法流程图因为其简单但又是求解其他方程的基础，所以发展了许多典型的解法.所有算法中的f 就是代表上式中(,)f x y ，而y f 表示(,)f x y y ∂∂，x f 表示(,)f x y x∂∂. 简单欧拉法是一种单步递推算法.简单欧拉法的公式如下所示：1(,)n n n n y y hf x y +=+ (21)简单欧拉法的算法过程介绍如下:给出自变量x 的定义域[,]a b ，初始值0y 及步长h .对0,1,()/k b a h =-，计算1(,)k k k k y y hf x y +=+5.2.2欧拉数值算法（隐式）隐式欧拉法也叫退欧拉法，隐式欧拉法的公式如下所示：111(,)n n n n y y hf x y +++=+ (22)隐式欧拉法是一阶精度的方法，比它精度高的公式是：111[(,)(,)]2n n n n n n hy y f x y f x y +++=++ (23)隐式欧拉的算法过程介绍如下.给出自变量x 的定义域[,]a b ，初始值0y 及步长h .对0,1,()/k b a h =-，用牛顿法或其他方法求解方111(,)k k k k y y hf x y +++=+得出1k y +.5.2.3 欧拉预估-校正法改进欧拉法是一种二阶显式求解法，其计算公式如下所示：1[,(,)]22n n n n n n h h y y hf x y f x y +=+++11(,)[(,)(,)]2n n n n n n n n t y hf x y h y y f x y f x t ++=+⎧⎪⎨=++⎪⎩ (24)四阶龙格-库塔法有多种形式，除了改进的欧拉法外还有中点法.中点法计算公式为：1[,(,)]22n n n n n n h h y y hf x y f x y +=+++ (25)5.3 模型求解5.3.1运用MATLAB 软件对模型求解结果及分析用欧拉法、改进的欧拉格式、经典的四阶龙格库塔法来求解常微分方程的边值问题，并且比较其精度（具体的MATLAB 源程序见附录） 以下进行实例分析：实例一. 1(0)1dyy x dx y ⎧=-++⎪⎨⎪=⎩由题可知精确解为：x y x e -=+，当x=0时，y(x)=0.在这里取步长h 为0.1， 通过MATLAB 程序的计算，相应的结果如下：表5-1 步长为0.1时各方法的预测值与精确值的比较(精确到6位小数)初值 Euler 法 相对误差 预估校正法 相对误差 经典四阶库 相对误差精确值 0 -- -- -- -- -- -- 1.00000 0.1 1.00910 0.00424 1.00500 0.00016 1.00484 0.00000 1.00484 0.2 1.02646 0.00759 1.01903 0.00029 1.01873 0.00000 1.01873 0.3 1.05134 0.01011 1.04122 0.00038 1.04082 0.00000 1.04082 0.4 1.08304 0.01189 1.07080 0.00045 1.07032 0.00000 1.07032 0.5 1.12095 0.01303 1.10708 0.00049 1.10653 0.00000 1.10653 0.6 1.16451 0.01366 1.14940 0.00052 1.14881 0.00000 1.14881 0.7 1.21319 0.01388 1.19721 0.00052 1.19659 0.00000 1.19659 0.8 1.26655 0.01378 1.24998 0.00052 1.24933 0.00000 1.24933 0.9 1.32414 0.01344 1.30723 0.00050 1.30657 0.00000 1.30657 1.01.38558 0.01294 1.36854 0.00048 1.36788 0.000001.36788步长为0.1时的精确值与预测值的比较精确值欧拉法改进欧拉格式四阶龙格库塔轴Y00.10.20.30.40.50.60.70.80.91 1.1 1.2 1.3 1.4X轴图5.2 步长为0.1时各方法的预测值与精确值的比较原函数图像轴YX轴图5.3 步长为0.1时原函数图像在这里取步长h为0.05，通过MATLAB程序的计算，相应的结果如下：表5-2 h=0.05时三个方法与精确值的真值表步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0 -- -- -- -- -- --1.00000 0.05 1.00250 0.00911 1.00125 0.01035 1.00123 0.01037 1.00484 0.10 1.00738 0.01711 1.00488 0.01954 1.00484 0.01958 1.01873 0.15 1.01451 0.02405 1.01076 0.02765 1.01071 0.02770 1.04082 0.20 1.02378 0.03001 1.01880 0.03473 1.01873 0.03479 1.07032 0.25 1.03509 0.03507 1.02889 0.04085 1.02880 0.04093 1.10653 0.30 1.04834 0.03930 1.04092 0.04610 1.04082 0.04619 1.14881 0.35 1.06342 0.04277 1.05480 0.05053 1.05469 0.05063 1.19659 0.40 1.08025 0.04555 1.07044 0.05422 1.07032 0.05432 1.24933 0.45 1.09874 0.04772 1.08775 0.05724 1.08763 0.05734 1.30657 0.50 1.11880 0.04933 1.10666 0.05964 1.10653 0.05975 1.36788 0.55 1.14036 0.05045 1.12709 0.06150 1.12695 0.06161 1.00484 0.60 1.16334 0.05113 1.14895 0.06286 1.14881 0.06298 1.01873 0.65 1.18768 0.05143 1.17219 0.06379 1.17205 0.06391 1.04082 0.70 1.21329 0.05139 1.19674 0.06433 1.19659 0.06445 1.07032 0.75 1.24013 0.05106 1.22252 0.06453 1.22237 0.06465 1.10653 0.80 1.26812 0.05048 1.24949 0.06443 1.24933 0.06455 1.14881 0.85 1.29722 0.04969 1.27757 0.06408 1.27742 0.06419 1.19659 0.90 1.32735 0.04871 1.30673 0.06349 1.30657 0.06361 1.249330.95 1.35849 0.04759 1.33690 0.06272 1.33674 0.06283 1.306571.00 1.39056 0.04634 1.36804 0.06178 1.36788 0.06189 1.36788欧拉格式 改进欧拉格式四阶龙格库塔精确值图5.4 步长为0.05时各方法的预测值与精确值的比较11.051.11.151.21.251.31.351.4X 轴Y 轴原函数图像图5.5 步长为0.05时原函数图像实例2 (0)1xdy xe ydx y -⎧=-⎪⎨⎪=⎩由题可知精确解为：21(2)2x x e -+当x=0时，y(x)=0.在这里取步长h为0.1，通过MATLAB程序的计算，相应的结果如下：步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0.1 0.9000 0.0318 0.9096 0.0214 0.9094 0.0216 0.9295 0.2 0.8192 0.0611 0.8359 0.0420 0.8356 0.0424 0.8726 0.3 0.7544 0.0876 0.7761 0.0614 0.7757 0.0619 0.8268 0.4 0.7027 0.1109 0.7277 0.0793 0.7272 0.0799 0.7903 0.5 0.6617 0.1310 0.6886 0.0956 0.6881 0.0963 0.7615 0.6 0.6294 0.1480 0.6572 0.1103 0.6567 0.1110 0.7387 0.7 0.6040 0.1621 0.6320 0.1234 0.6315 0.1241 0.7209 0.8 0.5841 0.1737 0.6116 0.1349 0.6111 0.1355 0.70690.9 0.5686 0.1830 0.5951 0.1450 0.5946 0.1456 0.69591.0 0.5563 0.1905 0.5815 0.1538 0.5811 0.1544 0.6872原函数图像轴Y00.10.20.30.40.50.60.70.80.91X轴图5.6 步长为0.1时原函数图像各方法的预测值与精确值的比较欧拉格式改进的格式四阶龙格库塔精确值轴Y0.10.20.30.40.50.60.70.80.91X轴图5.7 步长为0.1时各方法的预测值与精确值的比较在这里取步长h为0.05，通过MATLAB程序的计算，相应的结果如下：表5-4 步长为0.1时各方法的预测值与精确值的比较(精确到5位小数)步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0.050.95000 0.01342 0.95245 0.01088 0.95242 0.01090 0.96292 0.100.90490 0.02650 0.90947 0.02158 0.90942 0.02163 0.92953 0.150.86428 0.03916 0.87067 0.03206 0.87061 0.03213 0.89951 0.200.82774 0.05137 0.83567 0.04228 0.83559 0.04237 0.87256 0.250.79491 0.06307 0.80414 0.05219 0.80405 0.05230 0.84842 0.300.76545 0.07423 0.77576 0.06176 0.77566 0.06188 0.82682 0.350.73904 0.08482 0.75023 0.07096 0.75013 0.07110 0.80754 0.40 0.71541 0.09482 0.72730 0.07977 0.72719 0.07991 0.79035 0.45 0.69427 0.10422 0.70672 0.08817 0.70660 0.08832 0.77505 0.50 0.67539 0.11302 0.68825 0.09615 0.68813 0.09630 0.76146 0.55 0.65855 0.12123 0.67168 0.10370 0.67156 0.10386 0.74940 0.60 0.64352 0.12886 0.65683 0.11084 0.65671 0.11100 0.73871 0.65 0.63012 0.13593 0.64351 0.11756 0.64340 0.11773 0.72925 0.70 0.61819 0.14245 0.63157 0.12388 0.63145 0.12405 0.72087 0.75 0.60754 0.14847 0.62085 0.12981 0.62073 0.12998 0.71347 0.80 0.59805 0.15400 0.61121 0.13537 0.61110 0.13553 0.70691 0.85 0.58957 0.15908 0.60254 0.14058 0.60243 0.14073 0.70109 0.90 0.58198 0.16374 0.59470 0.14545 0.59460 0.14560 0.695930.95 0.57517 0.16802 0.58761 0.15001 0.58751 0.15016 0.691321.00 0.56904 0.17194 0.58117 0.15429 0.58107 0.15443 0.687190.10.20.30.40.50.60.70.80.910.550.60.650.70.750.80.850.90.951X 轴Y 轴各方法下的预测值与精确值的比较欧拉格式改进欧拉格式四阶龙格库塔精确值图5.8 步长为0.05时各方法的预测值与精确值的比较六、模型的评价本文着重讨论了4阶的隆格库塔法来求解微分方程，并且通过两个实例验证了隆格库塔法在求解初值问题的优越性.从上面的实例比较可知，在计算精度上，四阶经典龙格-库塔方法的误差最小，改进欧拉方法其次，欧拉方法误差则比较大，所以四阶经典龙格-库塔方法得到最佳的精度.而在计算量上面，相应地，很明显的四阶经典龙格-库塔方法也是最大，改进欧拉方法其次，欧拉方法计算量最小.这样的结果，说明了运用以上三种方法时，其计算量的多少与精度的大小成正比.我们在实际运用与操作中，可以根据实际情况，选择这3种方法中的其中一种最适合的，追求精度的话，可以使用四阶经典龙格-库塔方法；而改进的欧拉方法，在精度上和计算量上都表现得很出色，能够满足一般情况；而欧拉方法更主要的是适用于对y的估计上，而精度则有所欠缺，以上各方法的选择，都取决于具体的情况.七、课程设计的总结与体会本文着重采用隆格库塔法运用MATLAB编程来求解微分方程，相比于欧拉法以及预估校正法，隆格库塔法在提高近似值解的精度上是非常起作用的，而且又具有计算量不大、算法组织容易.其次，每一次的课程设计总是让我学到了更多的知识，不论是C++、SPASS还是MATLAB软件，这些都让我学到了如何解决实际问题的好工具，通过这些工具，是自己能够得到突破和成长.以下是我完成此次课程设计的几点体会：（1）必须学好基础知识，在做的过程中，发现自己有很多东西都不懂，要博学必须从一点一点做起.以往训练得少只是把握的不牢靠.所以做起来感到有点吃力.所以，无论什么学科，一定要打好基础.（2）程序设计要靠多练，多见识，那样形成一种编程思维，我想对我是有很大好处的.尤其像我这种平时学得不扎实的人.（3）做事情要有恒心，遇到困难不要怕，坚决去做.如果做出来了，固然高兴，如果没有做出来也没关系，自己努力了对得起自己就好.同时，把它看做是对自己的锻炼. （4）做程序特别是做大程序是很有趣的.虽然有的问题很难，要花很多时间很多精力，但是那种解决了一个问题时的喜悦足以把付出的辛苦补偿回来.得到一种心里的慰藉.参考文献[1] 李庆杨，王能超，易大义编．数值分析（第四版）[M]：华中科技大学出版社，2006.[2] 姜启源，谢金星，叶俊编．数学模型（第三版）[M]．北京：高等教育出版社，2005[3] 刘琼荪，数学实验[M]，高等教育出版社，2004[4] 王建伟，MATLAB7.X程序设计[M]，中国水利水电出版社，2007[5] 王高雄，周之铭等编.常微分方程（第三版） [M]:高等教育出版社，2006[6]何坚勇编著. 运筹学基础（第二版）[M]. 北京：清华大学出版社，2008.附录附录一：显示欧拉法matlab程序%欧拉法clear allclcx=[];y=[];y1=[];h=0.1;x=0:h:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold ony1(1)=1;for j=2:ny1(j)=y1(j-1)+h*f(x(j-1),y1(j-1)); endY=[x;y1];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);迭代n次的后退的欧拉格式matlab程序%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);%迭代n次的后退的欧拉格式clear allclcN=input('请输入迭代次数:');x=[];y=[];y1=[];h=0.1;x=0:0.1:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');for j=2:ny1(j)=y1(j-1)+h*f(x(j-1),y1(j-1)); T=y1(j);for k=1:NT=y1(j-1)+h*f(x(j),T);endy1(j)=T;endY=[x;y1];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y); fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)附录二：预估校正法matlab程序%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);%预估校正法%欧拉法clear allclcx=[];y1(1)=1;y2=[];y2(1)=1;h=0.1;x=0:0.1:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold onfor j=2:ny1(j)=y2(j-1)+h*f(x(j-1),y2(j-1));y2(j)=y2(j-1)+(h/2).*[f(x(j-1),y2(j-1))+f(x(j),y1(j))]; endY=[x;y2];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y2,'r-');figure(2)DT=y-y2;plot(x,DT)附录三：四阶龙格库塔法matlab程序%四阶龙格库塔法clear allclcn=length(x);y=[];y1=[];y1(1)=1;for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold onfor j=2:nK1=f(x(j-1),y1(j-1));K2=f(x(j-1)+h/2,y1(j-1)+h/2*K1);K3=f(x(j-1)+h/2,y1(j-1)+h/2*K2);K4=f(x(j-1)+h,y1(j-1)+h*K3);y1(j)=y1(j-1)+(h/6)*(K1+2*K2+2*K3+K4); endY=[x;y1];fid=fopen('data1.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)。
常微分方程数值解法-欧拉法、改进欧拉法与四阶龙格库塔法常微分方程数值解法
y( xn1)
y( xn
Байду номын сангаас
h)
y(xn )
hy'( xn )
h2 2!
y''( )
进一步: 令
h2 y( xn ) hy'( xn ) 2! y''( xn )
常微分方 yn1 y( xn1 ) , yn y( xn )
程数值解
法-欧拉法 yn1 yn hf ( xn , yn ) h2
、改进欧 y( xn1 ) yn1
2
max y''( x)
a xb
拉法和四
三、Euler方法
已 知 初 值 问 题 的 一 般 形式 为:
dy
dx
f (x, y)
a xb
(1)
y( x0 ) y0
常微分方 用差商近似导数 程数值解 问题转化为
yn1 yn dy
h
dx
法-欧拉法 yn1 yn hf ( xn , yn )
法-欧 y(拉0) 法1
、改进欧
拉法和四
四、几何意义
由 x0 , y0 出发取解曲线 y yx 的切线（存在！），则斜率
dy
f x0, y0
dx x y
,
0
0
常微分方 由于 f x0, y0 及 x0, y0 已知，必有切线方程。
由点斜式写出切程线方数程：值解
法、-改欧进拉欧法 ddxy y y0 x x0
常微分方 程数值解 能用解析方法求出精确解的微分方程为数不多,
而且有的方程即使有解析解,也可能由于解的表达
法-欧拉法 式非常复杂而不易计算,因此有必要研究微分方程
龙格 库塔方法 数学课的讲解
得到高精度方法的一个直接想法是利用Taylor展开
假设式
y' =f(x,y) (a≤x≤b)
中的 f(x,y) 充分光滑,将y(xi+1)在x i点作Taylor展开，若 取右端不同的有限项作为y(xi+1)的近似值，就可得到 计算y(xi+1)的各种不同截断误差的数值公式。
14
注意到，a2 b21 1, c1就是c2二阶12龙格 - 库塔公式，也就
是改进的欧拉法。
y
K2
yi 1
yi
1 2
K1
K2
K1 hf (xi , yi )
y y(x)
K1
K
K2 hf (xi h, yi K1)
xi
xi1
x
因此，凡满足条件式有一簇形如上式的计算格
式，这些格式统称为二阶龙格—库塔格式。因此改
3
Runge-Kutta 方法是一种高精度的单步法,简称R-K法 9.4.1 龙格-库塔(R-K)法的基本思想
Euler公式可改写成
yi1 yi K K hf (xi , yi )
则 yi+1 的 表 达 式 与 y(xi+1) 的 Taylor 展 开 式 的 前 两 项 完全相同,即局部截断误差为O(h2)。
答案是否定的！无论四个参数怎样选择,都不能使公式 的局部截断误差提高到三阶。
这说明每一步计算两个函数值的二阶R-K方法最高阶为 二阶。若要获得更高阶得数值方法,就必须增加计算函数 值的次数。
2020/3/22
17
9.4.3 三阶龙格—库塔法
为进一步提高精度，在区间[xi, xi+1]上除两点xi和
4阶经典龙格库塔公式求解微分方程
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。
四阶龙格库塔法(Runge-Kutta)求解微分方程
四阶龙格库塔法（Runge-Kutta ）求解微分方程张晓颖（天津大学 材料学院 学号：1012208027）1 引言计算传热学中通常需要求解常微分方程。
这类问题的简单形式如下：{),(')(00y x f y y x y == （1）虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中的多数微分方程需要采用数值解法求解。
初值问题（1）的数值解法有个基本特点，它们采取“步进式”，即求解过程顺着节点排序一步一步向前推进。
这类算法是要给出用已知信息y n 、 y n −i ……计算y n +1的递推公式即可。
2 龙格库塔法（Runge-Kutta ）介绍假设对于初值问题（1）有解 y = y (x ) ，用 Taylor 展开有：......)(!3)(!2)()()(321+'''+''+'+=+n n n n n x y h x y h x y h x y x y （2） 龙格库塔法（Runge-Kutta ）实质上是间接的使用 Taylor 级数法的一种方法。
对于差商hx y x y n n )()(1-+，根据微分中值定理，存在 0 < θ < 1 ，使得：)()()(1h x y hx y x y n n θ+'=-+ （3）于是对于 y = y (x ) ，可得到：))(,()()(1h x y h x hf x y x y n n n n θθ+++=+ （4）设))(,(*h x y h x f K n n θθ++=，为区间 [x n , x n +1 ] 上的平均斜率。
四阶龙格库塔格式中的*K 由下式计算得到：⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+),()2,2()2,2(),()22(6342312143211hK y h x f K K h y h x f K K h y h x f K y x f K K K K K h y y n n n n nn n n n n （5） 四阶龙格库塔法（Runge-Kutta ）的每一步需要四次计算函数值f ，其截断误差更低，计算的精度更高。
matlab_常微分方程数值解法
dt 2
简朴问题可以求得解析解，多数实际问题靠数值求解 。
第4页
一阶常微分方程(ODE )初值问题 ： ODE ：Ordinary Differential Equation
dy
f
(x,
y)
dx
x0 x xn
y(x0 ) y0
数值解法就是求y(x)在某些分立旳节点 xn 上旳近似值 yn，用以近似y(xn)
x0
y0
x1 f y(x), x dx
x0
x2 f y(x), x dx
x1
y(x1) f y(x1), x1 h
第17页
同样，在[x0,xn+1] ，积分采用矩形近似，得：
y(xn1) y0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
yn y(xn )
第5页
2、欧拉近似办法
2.1 简朴欧拉(L.Euler, 1707-1783)办法。
dy
dx
f
(y, x)
y(x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解
就是从初值开始，后一种函数值由前一种函数值得到。核 心是构造递推公式。
y0 y1 y2 yn
第6页
i 1,2,...
第36页
没有一种算法可以有效地解决所有旳 ODE 问题，因此 MATLAB 提供了多种ODE函数。
函数 ODE类
特点
阐明
型
ode45
非刚性 单步法；4，5 阶 R-K 措施；合计 大部分场合旳首选措施
截断误差为 (△x)3
ode23
非刚性 单步法；2，3 阶 R-K 措施；合计 使用于精度较低旳情形
四阶Runge-Kutta方法
2012-2013（1）专业课程实践论文四阶Runge-Kutta方法李元东，0818180107，R数学08-1班常微分方程初值问题的数值解法是求方程（1）的解在点列1(0,1,)n n n x x h n -=+= 上的近似值n y ，这里n h 是1n x -到n x 的步长，一般略去下标记为h 。
0(,)()dy f x y dxy x y ⎧=⎪⎨⎪=⎩ （1）经典的R K -方法是一个四阶的方法，它的计算公式是：112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ （2） R K -方法的优点是：单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值f 。
在用龙格库塔方法时,要注意n 的选择要合适，n 太大，会使计算量加大，n 太小，h 较大，可能会使误差增大。
因此选择合适的n 很重要。
我们要在考虑精度的基础上，选择合适的n 。
在此，用C 语言实现了龙格库塔方法。
四阶Runge-Kutta流程图#include<stdlib.h>#include<stdio.h>int RungeKutta(double y0,double a,double b,int n,double *x,double *y,double (*function)(double,double)){double h=(b-a)/n,k1,k2,k3,k4;int i;x[0]=a;y[0]=y0;for(i=0;i<n;i++){x[i+1]=x[i]+h;k1=function(x[i],y[i]);k2=function(x[i]+h/2,y[i]+h*k1/2);k3=function(x[i]+h/2,y[i]+h*k2/2);k4=function(x[i]+h,y[i]+h*k3);y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;}return 1;}double function(double x,double y){return y-2*x/y;}//例子求y'=y-2*x/y(0<x<1);y0=1;int main(){ int i;double x[6],y[6];printf("用四阶龙格-库塔方法\n");RungeKutta(1,1,2,5,x,y,function);for(i=0;i<6;i++)printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);return 1;}例1．求解方程⎩⎨⎧=<<+='110,0y x y x y 步长2.0=h ; 解：将程序中y x y return/*2-改为y x return + 1. 输入n b a y ,,,0的值：1；0；1；52. 显示输出结果：000000.00=x 时 ，000000.10=y ； 200000.01=x 时 ，242800.11=y ； 400000.02=x 时 ，583636.12=y ； 600000.03=x 时 ，044213.23=y ； 800000.04=x 时 ，651042.24=y ； 000000.15=x 时 ，436502.35=y ；例2．求解方程⎩⎨⎧=<<+=',1,10),1/(30y x x y y 步长2.0=h ； 解：将程序中y x y return /*2-改为()x y return +1/*31. 输入n b a y ,,,0的值：1；0；1；52. 显示输出结果：000000.00=x 时 ，000000.10=y ； 200000.01=x 时 ，727548.11=y ； 400000.02=x 时 ，742951.22=y ； 600000.03=x 时 ，094181.43=y ； 800000.04=x 时 ，829211.54=y ； 000000.15=x 时 ，996012.75=y ；。
MATLAB常微分方程数值解——欧拉法、改进的欧拉法与四阶龙格库塔方法
MATLAB常微分⽅程数值解——欧拉法、改进的欧拉法与四阶龙格库塔⽅法MATLAB常微分⽅程数值解作者：凯鲁嘎吉 - 博客园1.⼀阶常微分⽅程初值问题2.欧拉法3.改进的欧拉法4.四阶龙格库塔⽅法5.例题⽤欧拉法，改进的欧拉法及4阶经典Runge-Kutta⽅法在不同步长下计算初值问题。
步长分别为0.2,0.4,1.0.matlab程序：function z=f(x,y)z=-y*(1+x*y);function R_K(h)%欧拉法y=1;fprintf('欧拉法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K=f(x,y);y=y+h*K;fprintf('欧拉法：x=%f, y=%f\n',x+h,y);endfprintf('\n');%改进的欧拉法y=1;fprintf('改进的欧拉法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K1=f(x,y);K2=f(x+h,y+h*K1);y=y+(h/2)*(K1+K2);fprintf('改进的欧拉法：x=%f, y=%f\n',x+h,y);endfprintf('\n');%龙格库塔⽅法y=1;fprintf('龙格库塔法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K1=f(x,y);K2=f(x+h/2,y+(h/2)*K1);K3=f(x+h/2,y+(h/2)*K2);K4=f(x+h,y+h*K3);y=y+(h/6)*(K1+2*K2+2*K3+K4);fprintf('龙格库塔法：x=%f, y=%f\n',x+h,y);end结果：>> R_K(0.2)欧拉法：x=0.000000, y=1.000000欧拉法：x=0.200000, y=0.800000欧拉法：x=0.400000, y=0.614400欧拉法：x=0.600000, y=0.461321欧拉法：x=0.800000, y=0.343519欧拉法：x=1.000000, y=0.255934改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=0.200000, y=0.807200改进的欧拉法：x=0.400000, y=0.636118改进的欧拉法：x=0.600000, y=0.495044改进的欧拉法：x=0.800000, y=0.383419改进的欧拉法：x=1.000000, y=0.296974龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=0.200000, y=0.804636龙格库塔法：x=0.400000, y=0.631465龙格库塔法：x=0.600000, y=0.489198龙格库塔法：x=0.800000, y=0.377225龙格库塔法：x=1.000000, y=0.291009>> R_K(0.4)欧拉法：x=0.000000, y=1.000000欧拉法：x=0.400000, y=0.600000欧拉法：x=0.800000, y=0.302400改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=0.400000, y=0.651200改进的欧拉法：x=0.800000, y=0.405782龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=0.400000, y=0.631625龙格库塔法：x=0.800000, y=0.377556>> R_K(1)欧拉法：x=0.000000, y=1.000000欧拉法：x=1.000000, y=0.000000改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=1.000000, y=0.500000龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=1.000000, y=0.303395注意：在步长h为0.4时，要将for i=1:1/h改为for i=1:0.8/h。
龙格库塔方法
• 可以通过检查步长折半前后两次计
算结果的偏差
=
(h)
y2 n1
y(h) n1
变步长方
来判断所选取的步长是否合法适
（1）对于给定的精度，如果 ,则反复将步长折半 进行计算直至 为止，这时取折半以后的“ 新值” 作为结果；
四阶RungeKutta方法
这样，下一个值(yn+1)由现在的值(yn)加上时间间隔 (h)和一个估算的斜率的乘积 决定。该斜率是以下斜率的加权平均：
k1是时间段开始时的斜率； k2是时间段中点的斜率，通过欧拉法采用斜率 k1 来 决定 y在点 tn + h/2的值； k3也是中点的斜率，但是这次采用斜率 k2决定 y值； k4是时间段终点的斜率，其 y值用 k3 决定。 当四个斜率取平均时，中点的斜率有更大的权值：
75
0.0033
误差 0.45e-4 0.17e-4 0.15e-4 0.48e-4 0.25e-4 0.55e-4 0.14e-4
y(xn) 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492
改进Euler法一步需要计算两个函数值 (h=0.1) 四阶Runge-Kutta方法一步需要计算四 个函数值（h=0.2） 总计算量大致相当，但四阶RungeKutta方法精度更高
设从节点xn出发，先以h为步长求出一个近似值，记为yn(h)1 因为经典格式的局部截断误差为O(h5),因此有
y(xn+1)-y(nh)1 Ch5 其中C与y(5)(x)在[xn, xn1]内的值有关
将步长折半，即取
h 2
为步长从xn跨两步到xn+1，求得一个近似值yn( h21)
Runge-Kutta法
3½×RK· ¨
1.422
¸ ĽøEuler· ¨
1.42
1.418 0.2965 0.297 0.2975 0.298 0.2985 0.299 0.2995 0.3 0.3005 17
§ 7.2 Runge-Kutta法
龙格-库塔(Runge-Kutta)方法简称R-K法，是一种应用较广的 高精度的单步法。
所谓单步法就是在计算yi时只用到前一步信息yi-1的方法。
本节介绍R-K法的构造原理、常用公式。
1
一、Runge-Kutta方法的构造原理 对于常微分方程的初值问题
y f ( x , y ) y( a ) y0 a xb
另一方面
h2 y( xi 1 ) y( xi ) hy( xi ) y( xi ) o(h3 ) 2
9
式（6）的局部截断误差:
en (h) y ( xi 1 ) yi 1 1 h 1 (1 2 ) y( xi ) h2 22 y( xi ) o(h3 ) 2
1 2 1 1 2 2 2
有无穷多组解，从而可以得到许多具体的二阶R -K公式 如：
1 2 1，2 1
1 1 0， 2 1， 2 2
10
用类似的方法也可以构造三阶R-K公式
h yi yi 1 ( K1 4 K2 K3 ) 6 K1 f ( xi 1, yi 1 ) h h K 2 f ( xi 1 , yi 1 K1 ) 2 2 K3 f ( xi 1 h, yi 1 h(2 K2 K1 ))
因而方法(8)有4阶精度
12
例1. 使用高阶R-K方法计算初值问题
四阶龙格库塔法微分方程求解发散
四阶龙格库塔法微分方程求解发散下载提示：该文档是本店铺精心编制而成的，希望大家下载后，能够帮助大家解决实际问题。
文档下载后可定制修改，请根据实际需要进行调整和使用，谢谢!本店铺为大家提供各种类型的实用资料，如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等，想了解不同资料格式和写法，敬请关注!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!四阶龙格库塔法微分方程求解发散引言微分方程是数学中一种重要的工具，用于描述自然界中许多现象的变化规律。
matlab用四阶龙格库塔函数求解微分方程组
一、介绍Matlab作为一种强大的科学计算软件，提供了众多函数和工具来解决微分方程组。
其中，四阶龙格库塔函数是一种常用的数值方法，用于求解常微分方程组。
本文将介绍如何使用Matlab中的四阶龙格库塔函数来求解微分方程组，并对该方法的原理和实现进行详细说明。
二、四阶龙格库塔方法四阶龙格库塔方法是一种常用的数值方法，用于求解常微分方程组。
它是一种显式的Runge-Kutta方法，通过逐步逼近微分方程的解，在每一步使用多个中间值来计算下一步的解。
该方法通过四个中间值来计算下一步的状态，并且具有较高的精度和稳定性。
三、在Matlab中使用四阶龙格库塔方法求解微分方程组在Matlab中，可以使用ode45函数来调用四阶龙格库塔方法来解决微分方程组的问题。
ode45函数是Matlab提供的用于求解常微分方程组的函数，可以通过指定微分方程组以及初值条件来调用四阶龙格库塔方法来进行求解。
1. 定义微分方程组我们需要定义要求解的微分方程组。
可以使用Matlab中的匿名函数来定义微分方程组，例如：```matlabf = (t, y) [y(2); -sin(y(1))];```其中，f是一个匿名函数，用于表示微分方程组。
在这个例子中，微分方程组是y' = y2, y2' = -sin(y1)。
2. 指定初值条件和求解区间接下来，我们需要指定微分方程组的初值条件和求解区间。
初值条件可以通过指定一个初始时刻的状态向量来完成，例如：```matlabtspan = [0, 10];y0 = [0, 1];```其中，tspan表示求解区间，y0表示初值条件。
3. 调用ode45函数进行求解我们可以通过调用ode45函数来求解微分方程组的数值解。
具体的调用方式如下：```matlab[t, y] = ode45(f, tspan, y0);```其中，t和y分别表示求解的时间点和对应的状态值。
四、示例下面我们通过一个具体的例子来演示如何使用Matlab中的四阶龙格库塔方法来求解微分方程组。
四阶Runge-Kutta法求解常微分方程
【SOS】在matlab中四阶Runge-Kutta法求解常微分方程匚悬赏分：150 -解决时间：2008-9-16 10:19【SOS】在matlab中四阶Runge-Kutta法求解常微分方程dx1/dt=x2 x1(0)=1e-8dx2/dt=x3 x2(0)=0dx3/dt=x4 x3(0)=0dx4/dt=f(t,x1,x2,x3,x4) x4(0)=0求解区间［0,1e-6］，在matlab中用四阶Runge-Kutta法求解，怎么用从一阶循环到四阶？有没有类似的具体算例？最好能帮小弟写个子程序希望各路高人出手相助，关乎小弟能否明年一月份顺利毕业，在此先多谢各位高人了。
提问者：zhihaikou -三级最佳答案没试过matlab，算这玩意太慢了，有fortran版的要不，有兴趣的话可以参考一下。
代码:SUBROUTINE run ge_kutta()!关于Runge-Kutta方法，该方法是用来解形如y'=f(t,y)的常微分方程的!经典的4阶R-K方法的公式如下：! Yn+1 = Yn + h/6 * (K1+2K2+2K3+K4)!其中! K仁f(Tn,Yn)! K2=f(T n+h/2,Y n+h/2*K1)! K3=f(T n+h/2,Y n+h/2*K2)! K4=f(T n+h,Y n+h*K3)type( no de_fish),po in ter :: pi7HO!S%（『!）siuo!」inu% 厂u 人 NOa%（r!）siuau;nu%r u 人NOd%（『!）siuo!」inu% 厂 u 人t7HN%（『!）siuo!」inu% 厂 u 人£ON%（『!）siuo!」inu%厂u 人（（>HO!S%lueu;nu _U!UJ >（ NOa%lueu;nu _u!iu > （ NOd%lueu;nu _u!iu >（ >HN%lueu;nu _U!UJ > eON%lueu;nu _u!iu >WW : >d31S -一ii（W+£>k0N+Z>k0N+"）*（09dois —|j ） + u 人二厂u 人L+u 人直具 *-以吐u 人甲：ed3丄S-—一ii（e>i*de;s -p+u 人'dois —：M+iuaurK ）—：|j ）uo!Qun 厂 |E 」60屮!=厂》£>UL|+u 人 POM u人'u+l DOMfeJB來：kSd3丄S ------- i（2>1*（0S/deis _p ）+UA 10S/deis _p+juajjno _p ）uo!pun4_|ej6aju!=e>j◊u ⑵q ）+u 人［1 變歪 u 人'乙/i|+i pOMfeJB Ux 來：c sd3丄s ------------- i（L>1*（0S/deis _|j ）+UA 107/deis _p+juajjno _p ）uo!pun4_|ej6aju!=3>j以*⑵L|）+u 人 H 變歪U 人'乙/1|+1 POMHB 辽X 來：2 2 C □丄s --------- i（uA^uejjno-pJuoipunPie 」6。
常微分方程数值解法及应用
常微分方程数值解法及其应用——浙江师范大学 数理信息工程学院【摘要】：本文对常微分方程初值问题现有的数值解法进行了综述研究。
主要讨论了几种常用的数值解法：即欧拉法，改进欧拉法，龙格库塔方法，阿达姆斯外插公式与内插公式等。
文章最后结合常见数值解法，对较为典型的微分方程模型进行数值求解，探讨了上述数值算法在实际建模问题中的应用。
【关键词】：常微分方程；数值解法；模型引言在工程技术问题中，经常需要求解常微分方程的初值问题()()00,dyf x y dxy x y ⎧=⎪⎨⎪=⎩（1） 而关于常微分方程各种各样的解析方法，只能求解一些特殊类型的方程。
在大 多数情况下，对初值问题(1)，只能用数值法求解。
数值解法的基本思想是求初值问题(1)的解()y x 在一系列等距节点：0123n x x x x x <<<<<处的近似值：0123,,,,,n y y y y y 。
其中相邻两个节点间的距离1i i h x x +=-称为步长，即节点0 (0,1,2,)k x x kh k n =+= 。
一、单步法单步法是指这类方法在计算1n y +时，只用到前一步的值1,,n n n x x y +，然后逐步往下计算。
这个算法的代表是龙格——库塔算法，简称R-K 方法。
四阶显示Runge —Kutta 方法是求解普通常微分方程初值问题数值解法中的重要方法，而隐式Runge —Kutta 公式是求解刚性常微分方程初值问题的重要方法。
（一）Euler 方法由微分由微分方程的基本概念可知，初值问题(1)的解是在xoy 平面上的一条过点()00,x y 的积分曲线()y y x =，在该曲线上任一点(),x y 处的切线斜率等于函数(),f x y 的值。
按照按照数值解法的基本思想，取等距节点0 (0,1,2,)k x x kh k n =+= ，在点()00,x y 处，以()00,f x y 为斜率作直线()()0000,y y f x y x x =+-，与直线1x x =交于点()11,x y 。
四阶龙格库塔截断误差
四阶龙格库塔截断误差四阶龙格库塔法是一种常用的数值方法，用于求解常微分方程的初值问题。
它的一个重要指标是截断误差，本文将以四阶龙格库塔截断误差为标题，探讨该误差的相关内容。
我们需要了解什么是截断误差。
在数值计算中，由于计算机的存储和运算能力有限，无法进行无限精度的计算。
因此，我们需要使用数值方法来近似求解数学问题。
而截断误差指的是数值方法得到的近似解与精确解之间的差异。
四阶龙格库塔法是一种经典的数值方法，用于求解常微分方程的初值问题。
它的基本思想是通过多次逼近来提高数值解的精度。
具体来说，四阶龙格库塔法将每一步的计算分为多个阶段，通过计算各个阶段的加权平均值来得到最终的数值解。
在四阶龙格库塔法中，每一步的计算可以分为四个阶段：计算函数值、计算斜率、计算加权平均值和更新步长。
在计算函数值阶段，我们利用已知的初始值和步长来计算出函数在当前点的值。
在计算斜率阶段，我们利用函数值来计算出当前点的斜率。
在计算加权平均值阶段，我们利用斜率来计算出当前点的加权平均值。
最后，在更新步长阶段，我们利用加权平均值和当前点的斜率来更新步长的大小。
通过多次逼近的过程，四阶龙格库塔法可以得到比其他低阶方法更精确的数值解。
但是，即使使用四阶龙格库塔法，我们仍然无法得到完全精确的解，因为截断误差是无法避免的。
截断误差的大小取决于步长的大小和函数的性质。
通常情况下，我们可以通过减小步长来减小截断误差的大小。
为了更好地理解截断误差的概念，我们可以通过一个简单的例子来说明。
假设我们要求解常微分方程dy/dx = x，在初始条件y(0) = 0的情况下，使用四阶龙格库塔法进行数值求解。
假设步长为h，我们可以得到数值解y_n。
然后，我们可以根据精确解y(x) = x^2/2来计算截断误差e_n = |y_n - y(x_n)|。
通过计算不同步长下的截断误差，我们可以观察到截断误差随着步长的减小而减小的趋势。
这是因为较小的步长意味着我们在每一步中进行更多的计算，从而得到更精确的数值解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
四阶龙格-库塔(R-K)方法求常微分方程中南大学MATLAB程序设计实践材料科学与工程学院2013年3月26日一、编程实现“四阶龙格－库塔（R-K ）方法求常微分方程”，并举一例应用之。
【实例】采用龙格-库塔法求微分方程：⎩⎨⎧==+-=0, 0)(1'00x x y y y1、算法说明：在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为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 n2、流程图：2.1、四阶龙格－库塔（R-K ）方法流程图：2.2、实例求解流程图：输入求解的求出待求用MATLAB 自带结束用自编函数四阶开始 输入待求微分方程、确定求解范围内k = 取否 求解：⎪⎪⎪⎩⎪⎪⎪⎨⎧++=++=++==) ,()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求解并输出：)22(63211k k k hy y nn +++=+ 是结束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、程序说明：利用平面应力状态下斜截面应力的一般公σxσyσyσxσατyxτyxτxyτxyτ式，画出任意平面应力状态下的应力圆（Moore圆），求出相应平面应力状态下的主应力（1σ、3σ），并求出该应力状态下任意方位角α的斜截面上的应力σ、τ。
2、程序流程图：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：'); end hold off%----------此部分求该应力状态下的主应力--------------------------------------------Sigma_Max=Smax Sigma_Min=Smin Tau_Max=Tmax Tau_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：0Sigma_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函数解方程组得到两椭圆的交点即方程组的解。