经典四阶龙格库塔法解一阶微分方程组
MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程
姓名:樊元君学号:02 日期:一、实验目的掌握MATLAB语言、C/C++语言编写计算程序的方法、掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。
掌握使用MATLAB程序求解常微分方程问题的方法。
:二、实验内容1、分别写出改进欧拉法与四阶龙格-库塔求解的算法,编写程序上机调试出结果,要求所编程序适用于任何一阶常微分方程的数值解问题,即能解决这一类问题,而不是某一个问题。
实验中以下列数据验证程序的正确性。
求,步长h=。
*2、实验注意事项的精确解为,通过调整步长,观察结果的精度的变化^)三、程序流程图:●改进欧拉格式流程图:~|●四阶龙格库塔流程图:]四、源程序:●改进后欧拉格式程序源代码:function [] = GJOL(h,x0,y0,X,Y)format longh=input('h=');…x0=input('x0=');y0=input('y0=');disp('输入的范围是:');X=input('X=');Y=input('Y=');n=round((Y-X)/h);\i=1;x1=0;yp=0;yc=0;for i=1:1:nx1=x0+h;yp=y0+h*(-x0*(y0)^2);%yp=y0+h*(y0-2*x0/y0);%·yc=y0+h*(-x1*(yp)^2);%yc=y0+h*(yp-2*x1/yp);%y1=(yp+yc)/2;x0=x1;y0=y1;y=2/(1+x0^2);%y=sqrt(1+2*x0);%fprintf('结果=%.3f,%.8f,%.8f\n',x1,y1,y);:endend●四阶龙格库塔程序源代码:function [] = LGKT(h,x0,y0,X,Y)。
format longh=input('h=');x0=input('x0=');y0=input('y0=');disp('输入的范围是:');"X=input('X=');Y=input('Y=');n=round((Y-X)/h);i=1;x1=0;k1=0;k2=0;k3=0;k4=0;for i=1:1:n~x1=x0+h;k1=-x0*y0^2;%k1=y0-2*x0/y0;%k2=(-(x0+h/2)*(y0+h/2*k1)^2);%k2=(y0+h/2*k1)-2*(x0+h/2)/(y0+h/2*k1);% k3=(-(x0+h/2)*(y0+h/2*k2)^2);%k3=(y0+h/2*k2)-2*(x0+h/2)/(y0+h/2*k2);% k4=(-(x1)*(y0+h*k3)^2);%k4=(y0+h*k3)-2*(x1)/(y0+h*k3);%…y1=y0+h/6*(k1+2*k2+2*k3+k4);%y1=y0+h/6*(k1+2*k2+2*k3+k4);%x0=x1;y0=y1;y=2/(1+x0^2);%y=sqrt(1+2*x0);%fprintf('结果=%.3f,%.7f,%.7f\n',x1,y1,y);end·end*五、运行结果:改进欧拉格式结果:;}四阶龙格库塔结果:步长分别为:和时,不同结果显示验证了步长减少,对于精度的提高起到很大作用,有效数字位数明显增加。
龙格库塔法解微分方程组
龙格库塔法解微分方程组引言微分方程组是数学中经常遇到的问题,在物理、工程和自然科学中都有广泛的应用。
为了求解微分方程组,我们需要利用数值方法来逼近解析解。
本文将介绍一种常用的数值方法——龙格库塔法(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)为止。
微分方程组的龙格库塔法对于一阶微分方程组的初值问题,我们可以将其转化为向量形式。
c++经典四阶龙格库塔法解一阶微分方程组
首先,龙格库塔法(Runge-Kutta method)是一种常用的数值积分算法,它可以用来数值地求解一阶常微分方程(ODE)。
对于一维ODE,经典的四阶龙格库塔法(RK4)是最常用的方法之一。
下面是一个使用C++实现经典四阶龙格库塔法解一阶微分方程组的示例代码:cpp#include <iostream>#include <vector>// 定义微分方程组dy/dx = f(x, y)std::vector<double> f(double x, const std::vector<double>& y) {std::vector<double> dy(y.size());// 此处为具体的微分方程组表达式,请根据实际需求来定义dy[0] = x * y[0]; // 示例中假设有一个一维微分方程y' = x*yreturn dy;}// 经典四阶龙格库塔法void rk4(double x0, double y0, double h, int numSteps) {double x = x0;double y = y0;for (int i = 0; i < numSteps; ++i) {std::vector<double> k1 = f(x, std::vector<double>{y});std::vector<double> k2 = f(x + h / 2, std::vector<double>{y + h * k1[0] / 2});std::vector<double> k3 = f(x + h / 2, std::vector<double>{y + h * k2[0] / 2});std::vector<double> k4 = f(x + h, std::vector<double>{y + h * k3[0]});y += h * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6;x += h;std::cout << "x: " << x << ", y: " << y << std::endl;}}int main() {double x0 = 0.0; // 初值xdouble y0 = 1.0; // 初值ydouble h = 0.1; // 步长int numSteps = 10; // 迭代次数rk4(x0, y0, h, numSteps);return 0;}请注意,这只是一个简单的示例代码,用于演示如何使用经典四阶龙格库塔法求解一维微分方程组。
四阶龙格库塔法(runge-kutta)求振动方程
四阶龙格库塔法(runge-kutta)求振动方程四阶龙格库塔法可以用于求解振动方程。
振动方程通常是一个二阶微分方程,可以转化为一个一阶微分方程组来使用四阶龙格库塔法求解。
假设我们要求解的振动方程为:m * d^2x/dt^2 + c * dx/dt + k * x = 0其中,m是质量,c是阻尼系数,k是弹性系数,x是位移,t 是时间。
首先,将振动方程转化为一阶微分方程组。
引入新的变量v来表示速度,则有:dx/dt = vdv/dt = (-c * v - k * x) / m接下来,使用四阶龙格库塔法求解该一阶微分方程组。
首先将时间区间分割成n个小区间,假设每个小区间的宽度为h。
则利用四阶龙格库塔法的公式可以得到如下更新规则:k1 = h * f(t, x, v)k2 = h * f(t + h/2, x + k1/2, v + l1/2)k3 = h * f(t + h/2, x + k2/2, v + l2/2)k4 = h * f(t + h, x + k3, v + l3)其中,f(t, x, v)表示的是微分方程组的右侧。
根据上述更新规则,可以得到新的位移和速度:x_new = x + (k1 + 2k2 + 2k3 + k4) / 6v_new = v + (l1 + 2l2 + 2l3 + l4) / 6重复以上步骤,可以一步一步地求解出整个时间区间内的位移和速度。
需要注意的是,初始条件x(0)和v(0)需要提供。
总结起来,使用四阶龙格库塔法求解振动方程的步骤如下:1. 将二阶微分方程转化为一阶微分方程组。
2. 设置初始条件x(0)和v(0)。
3. 将时间区间分割成n个小区间,计算每个小区间的宽度h。
4. 根据上述更新规则,逐步求解出整个时间区间内的位移和速度。
希望以上解答对您有所帮助!。
四阶龙格库塔法解微分方程
四阶龙格库塔法解微分方程一、四阶龙格库塔法解一阶微分方程①一阶微分方程: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]。
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
四阶龙格库塔法求解微分方程组题目matlab 在数学和工程学中,微分方程组求解是一个常见的问题。
其中,四阶龙格库塔法是一种常用的求解微分方程组的数值方法。
在本文中,我们将介绍如何使用 MATLAB 实现四阶龙格库塔法求解微分方程组。
首先,我们需要定义微分方程组。
假设我们要求解以下微分方程组:y1' = -2*y1 + y2y2' = -y1 + 2*y2可以将这个微分方程组写成向量形式:Y' = f(Y)其中,Y = [y1; y2]f(Y) = [-2*y1 + y2; -y1 + 2*y2]接下来,我们需要编写四阶龙格库塔法的代码。
具体步骤如下:1. 定义初始值 Y0 和时间间隔 dt。
2. 编写 f(Y) 的函数代码。
3. 使用四阶龙格库塔法计算下一个时间步的 Y 值。
4. 重复步骤 3 直到计算到指定的时间点。
以下是完整的 MATLAB 代码:% 定义初始值 Y0 和时间间隔 dtt0 = 0;tf = 10;dt = 0.01;Y0 = [1; 0];% 编写 f(Y) 的函数代码function dy = f(Y)dy = [-2*Y(1) + Y(2); -Y(1) + 2*Y(2)];end% 使用四阶龙格库塔法计算下一个时间步的 Y 值function Y = rk4(Y, f, dt)k1 = dt*f(Y);k2 = dt*f(Y + 0.5*k1);k3 = dt*f(Y + 0.5*k2);k4 = dt*f(Y + k3);Y = Y + (k1 + 2*k2 + 2*k3 + k4)/6;end% 计算 Y 的值Y = Y0;for t = t0:dt:tfY = rk4(Y, @f, dt);disp(Y);end以上代码将输出从 t0 到 tf 时间段内每个时间步的 Y 值。
您可以根据需要修改微分方程组和时间间隔等参数。
4阶RungeKutta法求解一阶常微分方程
《MATLAB语言及应用》大作业姓名:学号:学院:班级:题目编号:2013年10月134阶Runge-Kutta法求解一阶常微分方程。
一、Runge-Kutta法的数学理论龙格-库塔(Runge-Kutta)方法就是一种在工程上应用广泛的高精度单步算法。
由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。
该算法就是构建在数学支持的基础之上的。
龙格库塔方法的理论基础来源于泰勒公式与使用斜率近似表达微分,它在积分区间多预计算出几个点的斜率,然后进行加权平均,用做下一点的依据,从而构造出了精度更高的数值积分计算方法。
如果预先求两个点的斜率就就是二阶龙格库塔法,如果预先取四个点就就是四阶龙格库塔法。
一阶常微分方程可以写作:y'=f(x,y),使用差分概念。
(Yn+1-Yn)/h= f(Xn,Yn)推出(近似等于,极限为Yn')Yn+1=Yn+h*f(Xn,Yn)另外根据微分中值定理,存在0<t<1,使得Yn+1=Yn+h*f(Xn+th,Y(Xn+th))这里K=f(Xn+th,Y(Xn+th))称为平均斜率,龙格库塔方法就就是求得K的一种算法。
利用这样的原理,经过复杂的数学推导(过于繁琐省略),可以得出截断误差为O(h^5)的四阶龙格库塔公式:K1=f(Xn,Yn);K2=f(Xn+h/2,Yn+(h/2)*K1);K3=f(Xn+h/2,Yn+(h/2)*K2);K4=f(Xn+h,Yn+h*K3);Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6)二、Runge-Kutta的算法与流程图在龙格-库塔法中,四阶龙格-库塔法的局部截断误差约为0(h5),被广泛应用于解微分方程的初值问题。
其算法公式为:[]()()()()112341213243/622,/2,/2*,/2,/2*,,n n n n n n n n n n y y h k k k k k f x y k f x h y h k k f x h y h k k f x h y hk +⎧⎫=++++⎪⎪=⎪⎪⎪⎪=++⎨⎬⎪⎪=++⎪⎪⎪⎪=++⎩⎭流程图:(1)、四阶龙格-库塔方法流程图:(2)、实例求解流程图:三、Runge-Kutta的Matlab实现function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次就是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数)n=floor((b-a)/h);%求步数x(1)=a;%时间起点y(:,1)=y0;%赋初值,可以就是向量,但就是要注意维数for ii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii));k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);k4=ufunc(x(ii)+h,y(:,ii)+h*k3);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;%按照龙格库塔方法进行数值求解end四、Runge-Kutta的算例实现例:求解常微分方程d(y)/d(x)=-2*y+2*x*x+2*x ,0≤x≤0、5 ,y(0)=1、编辑m文件Untitled3:fun=inline('-2*y+2*x*x+2*x');[x,y]=ode45(fun,[0,0、5],1)>> Untitled3x =0、01250、02500、03750、0625 0、0750 0、0875 0、1000 0、1125 0、1250 0、1375 0、1500 0、1625 0、1750 0、1875 0、2000 0、2125 0、2250 0、2375 0、2500 0、2625 0、2750 0、2875 0、3000 0、31250、33750、35000、36250、37500、38750、40000、41250、42500、43750、45000、46250、47500、48750、5000 y =1、00000、97550、95190、9073 0、8864 0、8663 0、8471 0、8287 0、8112 0、7944 0、7785 0、7633 0、7489 0、7353 0、7224 0、7103 0、6989 0、6883 0、6783 0、6690 0、6605 0、6526 0、6454 0、63880、6277 0、6231 0、6191 0、6157 0、6130 0、6109 0、6093 0、6084 0、6080 0、6083 0、6091 0、6104 0、6124 0、6148 0、6179。
四阶龙格——库塔法
四阶龙格——库塔法2013-2014(1)专业课程实践论文题目:四阶龙格—库塔法2'''()11()()()......()()2!!pp p n n n n n h h y y x hy x y x y x O h p ++=+++++ 进行离散化,所得计算公式必为p 阶方法,式中'''''()(,),()(,)(,)(,)....x y x f x y y x f x y f x y f x y ==++由此,我们能够想到,通过提高泰勒展开式的阶数,可以得到高精度的数值方法,从理论上讲,只要微分方程的解()y x 充分光滑,泰勒展开方法可以构造任意的有限阶的计算公式,但事实上,具体构造这种公式往往相当困难,因为符合函数(,())f x y x 的高阶导数常常是很烦琐的,因此,泰勒展开方法一般不直接使用,但是我们可以间接使用泰勒展开方法,求得高精度的计算方法。
首先,我们对欧拉公式和改进欧拉公式的形式作进一步的分析。
如果将欧拉公式和改进的欧拉公式改写成如下的形式:欧拉公式{111(,)n n n n y y hK K f x y +==+改进的欧拉公式11211()22n n y y h K K +=++, 1(,)n n K f x y =,21(,)n n K f x h y hK =++。
这两组公式都是用函数(,)f x y 在某些点上的值的线性组合来计算1()n y x +的近似值1n y +,欧拉公式每前进一步,就计算一次(,)f x y 的值。
另一方面它是1()n y x +在n x 处的一阶泰勒展开式,因而是一阶方法。
改进的欧拉公式每前进一步,需要计算两次(,)f x y 的值。
另一方面它在(,)n n x y 处的泰勒展开式与1()n y x +在n x 处的泰勒展开式的前三项完全相同,因而是二阶方法。
这启发我们考虑用函数(,)f x y 在若干点上的函数值的线性组合来构造计算公式。
经典四阶龙格库塔法解一阶微分方程组
1.经典四阶龙格库塔法解一阶微分方程组1.1运用四阶龙格库塔法解一阶微分方程组算法分析(1-1),(1-2)(1-3)(1-4)(1-5)(1-6)(1-7)(1-8)(1-9)(1-10)经过循环计算由推得……每个龙格-库塔方法都是由一个合适的泰勒方法推导而来,使得其最终全局误差为,一种折中方法是每次进行若干次函数求值,从而省去高阶导数计算。
4阶龙格-库塔方法(RK4)是最常用的,它适用于一般的应用,因为它非常精准,稳定,且易于编程。
1.2经典四阶龙格库塔法解一阶微分方程流程图图1-1 经典四阶龙格库塔法解一阶微分方程流程图1.3经典四阶龙格库塔法解一阶微分方程程序代码:#include <iostream>#include <iomanip>using namespace std;void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial[3], double resu[3],double h){double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1;t0=initial[0];x0=initial[1];y0=initial[2];f1=f(t0,x0,y0); g1=g(t0,x0,y0);f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2); g2=g(t0+h/2, x0+h*f1/2,y0+h*g1/2);f3=f(t0+h/2, x0+h*f2/2,y0+h*g2/2); g3=g(t0+h/2, x0+h*f2/2,y0+h*g2/2);f4=f(t0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3); x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6; resu[0]=t0+h;resu[1]=x1;resu[2]=y1;}int main(){double f(double t,double x, double y);double g(double t,double x, double y);double initial[3],resu[3];double a,b,H;double t,step;int i;cout<<"输入所求微分方程组的初值t0,x0,y0:";cin>>initial[0]>>initial[1]>>initial[2];cout<<"输入所求微分方程组的微分区间[a,b]:";cin>>a>>b;cout<<"输入所求微分方程组所分解子区间的个数step:";cin>>step;cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<<setprecision(10); H=(b-a)/step;cout<< initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl;for(i=0;i<step;i++){ RK4( f,g ,initial, resu,H);cout<<resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]<<endl;initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];}return(0);}double f(double t,double x, double y){double dx;dx=x+2*y;return(dx);}double g(double t,double x, double y){double dy;dy=3*x+2*y;return(dy);}1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示:应用所编写程序计算所给例题:其中初值为求解区间为[0,0.2]。
四阶龙格库塔法求解微分方程组
四阶龙格库塔法求解微分方程组微分方程是数学中的一个重要分支,它描述了自然界中许多现象的发展规律。
在实际应用中,我们经常需要求解微分方程组,以得到系统的演化规律和性质。
本文将介绍一种常用的求解微分方程组的数值方法——四阶龙格库塔法。
一、微分方程组的数值解法微分方程组是形如下式的方程集合: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$个点的自变量和因变量的值。
经典四阶龙格库塔法解一阶微分方程组
数值计算课程设计1. 经典四阶龙格库塔法解一阶微分方程组1.1 运用四阶龙格库塔法解一阶微分方程组算法分析x k 1 x k h(f 1 2 f 2 2 f 3 f 4)6hy k 1 y k h 6(g 1 2g 2 2g 3 g 4 )t k 1 t k h经过循环计算由 t 0,x 0, y 0推得 t 1,x 1,y1 t 2,x 2,y 2⋯⋯每个龙格 -库塔方法都是由一个合适的泰勒方法推导而来,使得其最终全局 误差为 O h N, 一种折中方法是每次进行若干次函数求值,从而省去高阶导数计 算。
4阶龙格-库塔方法 (RK4)是最常用的,它适用于一般的应用,因为它非常精 准,稳定,且易于编程。
f 1 f (t k , x k , y k ) ,hf 2 f (t k h 2,x khh2f 1,y k 2g 1)hf 3 f (t k h 2,x khf 2,y k 2g 2)f 4 f (t k h ,x k hf 3,y k hg 3)g 1 g(t k ,x k , y k ) h g 2g(t k 2,x khh2f 1,y k 2g 1)hg 3 g(t k 2, x khf 2h2,y k 2g 2) g 4 g(t k h,x k hf 3, y k hg 3)1-1)1-2)1-3)1-4)1-5)1-6)1-7)1-8)1-9)1-10 )经典四阶龙格库塔法解一阶微分方程组1.2 经典四阶龙格库塔法解一阶微分方程流程图1.3 经典四阶龙格库塔法解一阶微分方程程序代码:#include <iostream> #include <iomanip> using namespace std;void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial[3], double resu[3],double h){double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1; t0=initial[0];x0=initial[1];y0=initial[2]; f1=f(t0,x0,y0); g1=g(t0,x0,y0); f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2);g2=g(t0+h/2,x0+h*f1/2,y0+h*g1/2); f3=f(t0+h/2,x0+h*f2/2,y0+h*g2/2);g3=g(t0+h/2,x0+h*f2/2,y0+h*g2/2);f4=f(t0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3);图 1-1 经典四阶龙格库塔法解一阶微分方程流程图数值计算课程设计x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6;resu[0]=t0+h;resu[1]=x1;resu[2]=y1;}int main(){double f(double t,double x, double y);double g(double t,double x, double y);double initial[3],resu[3];double a,b,H;double t,step;int i;cout<<" 输入所求微分方程组的初值t0,x0,y0:"; cin>>initial[0]>>initial[1]>>initial[2]; cout<<" 输入所求微分方程组的微分区间[a,b]:"; cin>>a>>b;cout<<" 输入所求微分方程组所分解子区间的个数step:"; cin>>step;cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<<setprecision( 10);H=(b-a)/step;cout<< initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl;for(i=0;i<step;i++){ RK4( f,g ,initial, resu,H);cout<<resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]<<endl;initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];}return(0);}double f(double t,double x, double y){double dx;经典四阶龙格库塔法解一阶微分方程组dx=x+2*y; return(dx);}double g(double t,double x, double y) {double dy;dy=3*x+2*y; return(dy);}1.4 经典四阶龙格库塔法解一阶微分方程程序调试结果图示:图 1-2 经典四阶龙格库塔法解一阶微分方程算法程序调试图2. 高斯列主元法解线性方程组2.1 高斯列主元法解线性方程组算法分析使用伪代码编写高斯消元过程:for k=1 to n-1 dofor i=k+1 to n l<=a(i,k)/a(k,k) for j=k to n do a(i,j)<=a(i,j)-l*a(k,j) end %end of for j b(i)<=b(i)-l*b(k) end %end offor i end %end of for k最后得到 A ,b 可以构成上三角线性方程组 接着使用回代法求解上三角线性方程组应用所编写程序计算所给例题:其中初值为求解区间为 。
最新matlab 四阶龙格-库塔法求微分方程
m a t l a b四阶龙格-库塔法求微分方程Matlab 实现四阶龙格-库塔发求解微分方程从理论上讲,只要函数在某区间上充分光滑,那么它可以展开为泰勒级数,因此在该区间上的函数值可用各阶导数值近似地表示出来,反之其各阶导数值也可用某些函数值的线性组合近似地表示出来。
龙格-库塔法就是将待求函数)(t y 展开为泰勒级数,并用方程函数),(y f t 近似其各阶导数,从而迭代得到)(t y 的数值解。
具体来说,四阶龙格-库塔迭代公式为)22(6143211k k k k h n n ++++=+y y ),(1n n t k y f =)2/,2/(12hk h t k n n ++=y f)2/,2/(23hk h t k n n ++=y f),(33hk h t k n n ++=y f实验内容:已知二阶系统21x x= ,u x x x 5.02.04.0212+--= ,0)0()0(21==x x ,u 为单位阶跃信号。
用四阶龙格-库塔法求数值解。
分析步长对结果的影响。
实验总结:实验报告要求简要的说明实验原理;简明扼要地总结实验内容;编制m 文件,并给出运行结果。
报告格式请按实验报告模板编写。
进入matlab ,Step1:choose way1 or way2way1):可以选择直接加载M 文件(函数M 文件)。
way2):点击new ——function ,先将shier (函数1文本文件)复制运行;点击new——function,再将RK(函数2文本文件)运行;点击new——function,再将finiRK(函数3文本文件)运行;Step2:回到command页面输入下面四句。
[t,k]=finiRK45([0;0],150);%迭代150次,步长=20/150[t1 k1]=ode45(@shier,[0 -10],[0 0]);%调用matlab自带四阶龙格-库塔,对比结果[t2 k2]=ode45(@shier,[0 10],[0 0]);plot(t,k(1,:),'-',t1,k1(:,1),'*',t2,k2(:,1),'^')%在图形上表示出来补充:改变步长影响数据的准确性。
四阶龙格库塔法求解微分方程组
四阶龙格库塔法求解微分方程组四阶龙格-库塔法(RK4)是一种常用的数值求解微分方程组的方法。
它通过将微分方程转化为差分方程,通过逐步迭代来求得方程的近似解。
首先,考虑一个一阶微分方程组:dx/dt = f(x,t)其中,x是一个向量,f是一个向量函数,t是时间。
我们的目标是求解x(t)。
RK4法的基本思想是,在每个步骤中,通过预测下一个时间点的解并进行修正来逼近实际的解。
具体步骤如下:1.给定初始条件x(t0)和步长h,计算t=t0+h。
2. 计算k1 = hf(x(t0),t0),即利用初始点(x(t0),t0)计算出的斜率。
3. 计算k2 = hf(x(t0)+k1/2,t0+h/2),即利用中间点(x(t0)+k1/2,t0+h/2)计算出的斜率。
4. 计算k3 = hf(x(t0)+k2/2,t0+h/2),同理得到另一个中间点的斜率。
5. 计算k4 = hf(x(t0)+k3,t0+h),即利用最后一个中间点(x(t0)+k3,t0+h)计算的斜率。
6.根据k1,k2,k3和k4计算下一个点的值:x(t1)=x(t0)+(k1+2k2+2k3+k4)/67.重复步骤2到6,直到达到所需的终止时间。
接下来,让我们通过一个具体的例子来说明RK4方法的应用。
假设我们要求解如下的一阶微分方程组:dx/dt = f(x,y,z)dy/dt = g(x,y,z)dz/dt = h(x,y,z)其中f,g和h是已知的向量函数,我们需要找到x(t),y(t)和z(t)。
首先,我们需要将该微分方程组转化为差分方程组。
我们将离散时间点t0,t1,t2,...,tn代入微分方程,得到:x(t0+h)=x(t0)+h*f(x(t0),y(t0),z(t0))y(t0+h)=y(t0)+h*g(x(t0),y(t0),z(t0))z(t0+h)=z(t0)+h*h(x(t0),y(t0),z(t0))然后,我们通过逐步迭代的方式来求解方程组。
四阶龙格库塔法原理
四阶龙格库塔法原理四阶龙格库塔法是一种常用的数值解微分方程的方法,它在工程、物理、计算机科学等领域都有着广泛的应用。
四阶龙格库塔法的原理相对简单,但能够较为准确地求解微分方程,因此受到了广泛的关注和应用。
首先,我们来了解一下什么是微分方程。
微分方程是描述自然界中变化规律的数学模型,它包含了未知函数及其导数的方程。
在实际问题中,我们常常需要求解微分方程来得到系统的演化规律。
然而,有些微分方程并没有解析解,这时就需要借助数值方法来求解。
四阶龙格库塔法就是其中一种常用的数值方法。
四阶龙格库塔法的原理可以简单概括为以下几个步骤:首先,我们需要将微分方程化为一阶方程组的形式。
这是因为四阶龙格库塔法是针对一阶方程组的数值解法。
其次,我们需要选择一个合适的步长h。
步长的选择对于数值解的精度有着重要的影响,通常需要通过实验来确定一个合适的步长。
然后,我们可以通过迭代的方式来逐步求解微分方程。
四阶龙格库塔法通过多次迭代,利用当前点的斜率来估计下一个点的值,从而逐步逼近真实解。
最后,我们可以通过迭代得到的数值解来近似地描述微分方程的解。
四阶龙格库塔法在一定条件下能够较为准确地求解微分方程,因此在实际应用中有着广泛的价值。
四阶龙格库塔法的原理相对简单,但是需要注意的是,在具体应用时需要考虑步长的选择、迭代次数的确定等参数,以及对于不同类型的微分方程可能需要进行适当的调整。
此外,四阶龙格库塔法也有其局限性,对于某些特殊类型的微分方程可能并不适用,因此在具体应用时需要综合考虑。
总的来说,四阶龙格库塔法是一种常用的数值解微分方程的方法,它通过迭代的方式逐步逼近真实解,具有一定的精度和稳定性。
在实际应用中,我们需要根据具体问题的特点来选择合适的数值方法,并注意参数的选择和调整,以获得准确的数值解。
四阶龙格库塔法在工程、物理、计算机科学等领域都有着广泛的应用,对于研究和解决实际问题具有重要的意义。
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中的四阶龙格库塔方法来求解微分方程组。
龙格库塔方法求解微分方程例子
龙格库塔方法求解微分方程例子耦合非线性谐振子网络模型为研究物理和生物等众多领域提供了一个有用的工具。
特别的,在这个网络模型上会出现混沌现象。
而对于如何控制这种非线性引起的混沌现象以及如何使得这些谐振子达到同步化这些方面的研究,最近引起了极大的关注。
这是因为在现实生活中,我们所面临的大部分系统都是混沌的。
而混沌由于对初始条件具有极高的敏感性,使得我们无法预测这类系统(比如天气预报)。
所以如何控制混沌对于改进我们的生活有着很大的帮助。
本文主要讨论如何用无序来抑制混沌。
无序在物理系统中通常意味着破坏空间和时间的规则性。
然而近来研究表明,在非线性系统中,无序起着相反的作用,例如热涨落或者随机分布,它们可以导致某些规则的图案!我们考虑一一维耦合非线性阻力摆链,加入无序相位,其动力学方程如下:将动力学方程变为一阶微分方程组可得:其中为无序相位,均匀无规地从区间内取值。
我们采用自由边界条件:。
并取适当的参数:定义平均相速度:。
其中。
可以用来描述系统的宏观集体行为。
在数值模拟中,我们取,即求得20个时刻的平均速度。
接下来我们的任务是利用龙格库塔方法来数值求解微分方程组,并求得相应的平均相速度与无序范围k 之间的关系。
结果表明,当k 充分大时,不同时刻的平均相速度变得一致,即无序在某种程度上抑制了混沌!首先,我们来介绍龙格库塔方法。
考虑一函数,将在处泰勒展开,得到2'11sin sin()(2) (n=1,2...N)n n n n n n n ml mgl t q gq q t t w j k q q q ••+-+=-+++++-'11sin sin()(2)n nn n n n n n n t q f f gq q t t w j k q q q ••+-==--+++++-n j [,]k k p p -011011,,,N N N N q q q q f f f f ++===='1,1,0.75,0.7155,0.4,0.25,0.5l g g t t w k =======11()()Nn n jT jT Ns f ==å2/T p w =()jT s 60,61,...80j =()y t ()y t t +t[1]其中。
四阶经典龙格库塔例题
四阶经典龙格库塔例题经典的四阶龙格-库塔方法是一种数值积分算法,用于求解常微分方程的数值解。
下面是一个例题,我们将使用四阶龙格-库塔方法来求解:考虑如下的一阶常微分方程:dy/dx = x^2 + y^2,其中y(0) = 1。
我们的目标是求解在区间[0, 1]上的数值解。
为了使用四阶龙格-库塔方法,我们需要将问题转化为一个离散的迭代过程。
首先,我们将区间[0, 1]等分为n个小区间,每个小区间的宽度为h = (1-0)/n。
我们用xi表示第i个小区间的起始点,yi表示在xi处的数值解的近似值。
四阶龙格-库塔方法的迭代过程如下:1. 设置初始条件,x0 = 0,y0 = 1。
2. 对于i = 0, 1, 2, ..., n-1,执行以下步骤:a. 计算k1 = h (xi^2 + yi^2)。
b. 计算k2 = h ((xi + h/2)^2 + (yi + k1/2)^2)。
c. 计算k3 = h ((xi + h/2)^2 + (yi + k2/2)^2)。
d. 计算k4 = h ((xi + h)^2 + (yi + k3)^2)。
e. 计算yi+1 = yi + (k1 + 2k2 + 2k3 + k4)/6。
f. 更新xi+1 = xi + h。
3. 输出近似解y(x) = yn,其中n为等分的小区间数。
通过以上迭代过程,我们可以得到在区间[0, 1]上的数值解。
需要注意的是,选择合适的等分数n对于求解的精度很重要。
通常情况下,增加n的值可以提高数值解的精确度,但也会增加计算的复杂度。
希望以上解答能够满足你的要求。
如果你还有其他问题,欢迎继续提问。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
1.经典四阶龙格库塔法解一阶微分方程组1.1运用四阶龙格库塔法解一阶微分方程组算法分析(1-1),(1-2)(1-3)(1-4)(1-5)(1-6)(1-7)(1-8)(1-9)(1-10)经过循环计算由推得……每个龙格-库塔方法都是由一个合适的泰勒方法推导而来,使得其最终全局误差为,一种折中方法是每次进行若干次函数求值,从而省去高阶导数计算。
4阶龙格-库塔方法(RK4)是最常用的,它适用于一般的应用,因为它非常精准,稳定,且易于编程。
1.2经典四阶龙格库塔法解一阶微分方程流程图图1-1 经典四阶龙格库塔法解一阶微分方程流程图1.3经典四阶龙格库塔法解一阶微分方程程序代码:#include <iostream>#include <iomanip>using namespace std;void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial[3], double resu[3],double h){double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1;t0=initial[0];x0=initial[1];y0=initial[2];f1=f(t0,x0,y0); g1=g(t0,x0,y0);f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2); g2=g(t0+h/2, x0+h*f1/2,y0+h*g1/2);f3=f(t0+h/2, x0+h*f2/2,y0+h*g2/2); g3=g(t0+h/2, x0+h*f2/2,y0+h*g2/2);f4=f(t0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3); x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6; resu[0]=t0+h;resu[1]=x1;resu[2]=y1;}int main(){double f(double t,double x, double y);double g(double t,double x, double y);double initial[3],resu[3];double a,b,H;double t,step;int i;cout<<"输入所求微分方程组的初值t0,x0,y0:";cin>>initial[0]>>initial[1]>>initial[2];cout<<"输入所求微分方程组的微分区间[a,b]:";cin>>a>>b;cout<<"输入所求微分方程组所分解子区间的个数step:";cin>>step;cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<<setprecision(10); H=(b-a)/step;cout<< initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl;for(i=0;i<step;i++){ RK4( f,g ,initial, resu,H);cout<<resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]<<endl;initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];}return(0);}double f(double t,double x, double y){double dx;dx=x+2*y;return(dx);}double g(double t,double x, double y){double dy;dy=3*x+2*y;return(dy);}1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示:应用所编写程序计算所给例题:其中初值为求解区间为[0,0.2]。
计算结果为:图1-2 经典四阶龙格库塔法解一阶微分方程算法程序调试图2.高斯列主元法解线性方程组2.1高斯列主元法解线性方程组算法分析使用伪代码编写高斯消元过程:for k=1 to n-1 dofor i=k+1 to nl<=a(i,k)/a(k,k)for j=k to n doa(i,j)<=a(i,j)-l*a(k,j)end %end of for jb(i)<=b(i)-l*b(k)end %end of for iend %end of for k最后得到A,b可以构成上三角线性方程组接着使用回代法求解上三角线性方程组因为高斯消元要求a(k,k)≠0(k=1,2,3……n-1)这就需要对高斯消元过程进行完善,即使用高斯列主元法:其步骤为:①找主元:计算,并记录其所在行r ,②交换第r行与第k行;③以第k行为工具行处理以下各行,使得从第k列的第k+1行到第n行的元素全部为0;④得到增广矩阵的上三角线性方程组;⑤使用回代法对上三角线性方程组进行求解2.2高斯列主元法解线性方程组流程图图2-1 高斯列主元法解线性方程组流程图2.3高斯列主元法解线性方程组程序代码#include<iostream>#include <cmath>#define N 3using namespace std;void main(){int i,j,k,n,p;float t,s,m,a[N][N],b[N],x[N];cout<<"请输入方程组的系数"<<endl;for(i=0;i<N;i++){for(j=0;j<N;j++)cin>>a[i][j];}cout<<"请输入方程组右端的常数项:"<<endl;for(i=0;i<N;i++)cin>>b[i];for(j=0;j<N-1;j++){ p=j;for(i=j+1;i<N;i++) //寻找系数矩阵每列的最大值{if(fabs(a[i][j])>fabs(a[p][j]))p=i;}if(p!=j) //交换第p行与第j行{for(k=0;k<N;k++){t=a[j][k];a[j][k]=a[p][k];a[p][k]=t;} //交换常数项的第p行与第j行t=b[p];b[p]=b[j];b[j]=t;} //把系数矩阵第j列a[j][j]下面的元素变为0for(i=j+1;i<N;i++){ m=-a[i][j]/a[j][j];for(n=j;n<N;n++)a[i][n]=a[i][n]+a[j][n]*m;b[i]=b[i]+b[j]*m;}} //求方程组的一个解x[N-1]=b[N-1]/a[N-1][N-1]; //回代法求方程组其他解for(i=N-2;i>=0;i--){s=0.0;for(j=N-1;j>i;j--){s=s+a[i][j]*x[j];x[i]=(b[i]-s)/a[i][i];}}cout<<"方程组的解如下:"<<endl;for(i=0;i<=N-1;i++)cout<<x[i]<<endl;}2.4高斯列主元法解线性方程组程序调试结果图示:求解下列方程组图2-2 高斯列主元法解线性方程组程序算法调试图3.牛顿法解非线性方程组3.1牛顿法解非线性方程组算法说明牛顿法主要思想是用进行迭代。
因此首先需要算出的雅可比矩阵,再求过求出它的逆,当它达到某个精度时即停止迭代。
具体算法如下:首先设已知。
①:计算函数(3-1)②:计算雅可比矩阵(3-2)(3-3)③:求线性方程组的解。
④:计算下一点重复上述过程。
3.2牛顿法解非线性方程组流程图图3-1 牛顿法解非线性方程组流程图3.3牛顿法解非线性方程组程序代码#include<iostream>#include<cmath>#define N 2 // 非线性方程组中方程个数、未知量个数#define Epsilon 0.0001 // 差向量1范数的上限#define Max 100 //最大迭代次数using namespace std;const int N2=2*N;int main(){void ff(float xx[N],float yy[N]);//计算向量函数的因变量向量yy[N]void ffjacobian(float xx[N],float yy[N][N]);//计算雅克比矩阵yy[N][N] void inv_jacobian(float yy[N][N],float inv[N][N]);//计算雅克比矩阵的逆矩阵invvoid newdundiedai(float x0[N], float inv[N][N],float y0[N],float x1[N]);//由近似解向量 x0 计算近似解向量 x1floatx0[N]={2.0,0.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errorno rm;int i,j,iter=0;//如果取消对x0的初始化,撤销下面两行的注释符,就可以由键盘向x0读入初始近似解向量//for( i=0;i<N;i++)// cin>>x0[i];cout<<"初始近似解向量:"<<endl;for (i=0;i<N;i++)cout<<x0[i]<<" ";cout<<endl;cout<<endl;do{iter=iter+1;cout<<"第 "<<iter<<" 次迭代开始"<<endl;//计算向量函数的因变量向量 y0ff(x0,y0);//计算雅克比矩阵 jacobianffjacobian(x0,jacobian);//计算雅克比矩阵的逆矩阵 invjacobianinv_jacobian(jacobian,invjacobian);//由近似解向量 x0 计算近似解向量 x1newdundiedai(x0, invjacobian,y0,x1);//计算差向量的1范数errornormerrornorm=0;for (i=0;i<N;i++)errornorm=errornorm+fabs(x1[i]-x0[i]);if (errornorm<Epsilon) break;for (i=0;i<N;i++)x0[i]=x1[i];} while (iter<Max);return 0;}void ff(float xx[N],float yy[N]){float x,y;int i;x=xx[0];y=xx[1];yy[0]=x*x-2*x-y+0.5;yy[1]=x*x+4*y*y-4;cout<<"向量函数的因变量向量是: "<<endl;for( i=0;i<N;i++)cout<<yy[i]<<" ";cout<<endl;cout<<endl;}void ffjacobian(float xx[N],float yy[N][N]){float x,y;int i,j;x=xx[0];y=xx[1];//jacobian have n*n elementyy[0][0]=2*x-2;yy[0][1]=-1;yy[1][0]=2*x;yy[1][1]=8*y;cout<<"雅克比矩阵是: "<<endl;for( i=0;i<N;i++){for(j=0;j<N;j++)cout<<yy[i][j]<<" ";cout<<endl;}cout<<endl;}void inv_jacobian(float yy[N][N],float inv[N][N]) {float aug[N][N2],L;int i,j,k;cout<<"开始计算雅克比矩阵的逆矩阵:"<<endl;for (i=0;i<N;i++){ for(j=0;j<N;j++)aug[i][j]=yy[i][j];for(j=N;j<N2;j++)if(j==i+N) aug[i][j]=1;else aug[i][j]=0;}for (i=0;i<N;i++){ for(j=0;j<N2;j++)cout<<aug[i][j]<<" ";cout<<endl;}cout<<endl;for (i=0;i<N;i++){for (k=i+1;k<N;k++){L=-aug[k][i]/aug[i][i];for(j=i;j<N2;j++)aug[k][j]=aug[k][j]+L*aug[i][j];}}for (i=0;i<N;i++){ for(j=0;j<N2;j++)cout<<aug[i][j]<<" ";cout<<endl;}cout<<endl;for (i=N-1;i>0;i--){for (k=i-1;k>=0;k--){L=-aug[k][i]/aug[i][i];for(j=N2-1;j>=0;j--)aug[k][j]=aug[k][j]+L*aug[i][j];}}for (i=0;i<N;i++){ for(j=0;j<N2;j++)cout<<aug[i][j]<<" ";cout<<endl;}cout<<endl;for (i=N-1;i>=0;i--)for(j=N2-1;j>=0;j--)aug[i][j]=aug[i][j]/aug[i][i];for (i=0;i<N;i++){ for(j=0;j<N2;j++)cout<<aug[i][j]<<" ";cout<<endl;for(j=N;j<N2;j++)inv[i][j-N]=aug[i][j];}cout<<endl;cout<<"雅克比矩阵的逆矩阵: "<<endl;for (i=0;i<N;i++){ for(j=0;j<N;j++)cout<<inv[i][j]<<" ";cout<<endl;}cout<<endl;}void newdundiedai(float x0[N], float inv[N][N],float y0[N],float x1[N]) {int i,j;float sum=0;for(i=0;i<N;i++){ sum=0;for(j=0;j<N;j++)sum=sum+inv[i][j]*y0[j];x1[i]=x0[i]-sum;}cout<<"近似解向量:"<<endl;for (i=0;i<N;i++)cout<<x1[i]<<" ";cout<<endl;cout<<endl;}3.4牛顿法解非线性方程组程序调试结果图示:图3-2 牛顿法解非线性方程组程序算法调试图图3-3 牛顿法解非线性方程组程序算法调试图图3-4 牛顿法解非线性方程组程序算法调试图4.龙贝格求积分算法4.1龙贝格求积分算法分析生成j>=k 的近似积分结果逼近表,并以R(j+1,j+1)为最终解来逼近积分。