三九宝宝网宝宝教育学龄段教育

matlab实现欧拉法和RK4方法的数值计算

02月10日 编辑 39baobao.com

[如何实现合理的预算管控与费用报销]企业对于如何通过更有效的手段来实现合理的预算管控与费用报销,有着非常迫切的渴求。如何让企业的预算管控与费用报销不再是纸上谈兵,能有效配合企业的发展等一些列问题长期以...+阅读

程序已经写了,不过步长你得自己调,当步长较小时,计算时间会很长

另外,tend是时间的终值,你可以设小一些。因为解析解为10*cos(x),我设成pi,就是计算半个周期。

x''(t)=-x(t)

引入y1=x,y2=x',则

y1'=y2

y2'=-x=-y1

初始条件为:

y1(0)=10;

y2(0)=0;

将下面两行百分号之间的内容,保存成DiffEulerRk4.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function MaxDiffX=DiffEulerRk4(dt,PlotFlag)

%dt是步长

%PlotFlag是否作图

if nargin<1

dt=0.01;

end

if nargin<2

PlotFlag=0;

end

f=inline('[y(2);-y(1)]','t','y'); %微分方程的右边项

t0=0; %初始时刻

tend=pi; %计算的点数

tt=t0:dt:tend; %一系列离散的点

N=length(tt); %点数

y0=[10;0];

%%(1)欧拉法

EulerY=y0;

for i=2:N

EulerY(:,i)=EulerY(:,i-1)+dt*f(tt(i-1),EulerY(:,i-1));

end

EulerX=EulerY(1,:); %取x

%%(2)龙格库塔法

RkY=y0;

for i=2:N

k1=f(tt(i-1), RkY(:,i-1));

k2=f(tt(i-1)+dt/2, RkY(:,i-1)+k1*dt/2);

k3=f(tt(i-1)+dt/2, RkY(:,i-1)+k2*dt/2);

k4=f(tt(i-1)+dt, RkY(:,i-1)+k3*dt);

RkY(:,i)=RkY(:,i-1)+(k1+2*k2+2*k3+k4)*dt/6;

end

RkX=RkY(1,:); %取x

%精确解

syms t

analytic=dsolve('D2x=-x','x(0)=10','Dx(0)=0','t');

rightdata=subs(analytic,t,tt);

if PlotFlag

plot(tt,EulerX,'b-',tt,RkX,'r--',tt,rightdata,'g-.')

legend('Euler','Runge-Kutta','analytic')

end

MaxDiffX=[max(abs(RkX-rightdata)),max(abs(EulerX-rightdata))];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

所有题,都得你自己调步长。

输入:

DiffEulerRk4(0.01,1) %步长取0.01的计算结果,参数为1代表作图,自己得修改步长

%%下面是变化

Error=[];

Dt=[5e-4,1e-3,2e-3,5e-3,0.01,0.05,0.1];

for dt=Dt %几种步长,自行修改

dt %查看dt,步长小,计算量大

Error_1=DiffEulerRk4(dt); %不作图

Error=[Error;Error_1]; %保存欧拉法误差

end

semilogx(Dt,Error)

legend('Euler','RK4')

xlabel('步长')

ylabel('误差')

title('与理论值误差')

以下为关联文档:

LaTeX实现中文显示高手指导!比较齐全的宏包,复制过去就能用。 \documentclass[a4paper,4pt]{article} \usepackage{geometry}%页边距 \geometry{left=2.5cm,right=2.5cm,top=2.5cm,bottom=2.5cm} \lines...

如何在中学美术课堂教学中实现美术教育目标美育被明确地列入教育大纲,标志着我国社会主义教育思想的重大发展和社会主义教育制度的进一步完善。美术课的目的不是把所有学生培养成艺术创作家,但是对于提高学生的审美文化...

perl程序如何实现将一段文本插入某一行下方#!/usr/bin/perl#t.plopen A,"txt1"; #一段文本my a1=();while(){ #插入文本if (/I(19|12)/){print $_,"a1";}else {print $_}}===============perl t.pl txt2I19 (AGND AVDD vb_4...

perl解析报文每3行输出到1行怎么实现例如 a1 a2 a3 b1 b2 b3 c1用sed cat test_baidu a1 a2 a3 b1 b2 b3 c1 c2 c3 sed -e 'N;N;s/\n//g' test_baidu 用perl (perl 水平很烂 现学的) open('WCD','my $num = 1; my $str = ''; while() { c...

perl提取特定文本行的下一行如何实现open F, "F:/1.txt"; #1.txt是你画面上那个文件 open OUT, ">F:/out.txt"; #o.txt是你要保存的文件 array=; $count=-1; foreach (array){ $count++; if(/固定文本/){$start=$coun...

如何实现网络数据集的网络分析网络分析是arcgis提供的重要的空间分析的功能,利用它可以模拟现实世界的网络问题。例如多个地点的最短路径问题。 在arcgis中,将地理网络模型分为两种:几何网络模型和网络数据...

arcgis网络数据集最短路径是需要编程实现的么最短路径分析属于ArcGIS的 网络 分析范畴。而ArcGIS的网络分析分为两类,分别是基于几何网络和网络数据集的网络分析。它们都可以实现最短路径功能。下面先介绍基于几何网络的...

如何用matlab实现无限大功率电源供电系统的三相短路分析暂态上机计算教学指导书(适用于Matlab) 一、基础知识1. 为完成本次设计,需要掌握的Matlab的基础知识有:数组的创建和计算,循环语句,条件选择语句,条件判断语句,绘制二维图的语句;2....

页面调度算法的实现和分析源码dev c++ #include#include#includetypedef struct mem { int num; int v; }meme; static int process[1024],L,M; void LRU(); void FIFO(); void get(); int menu(); int m...

推荐阅读
图文推荐