matlab 欧拉法微分方程 MATLAB实战——微分方程组的解法之欧拉法与4阶龙格库塔法...

MATLAB实战——微分方程组的解法之欧拉法与4阶龙格库塔法

MATLAB实战——微分方程组的解法之欧拉法与4阶龙格库塔法

实现：

%%Question 1 part(b)

clear all;

clc;

t0 = 0;

x0 = 1/2;

dt = 0.1;

tf = 1;

t_range = t0:dt:tf;

x_EU = zeros(1,length(t_range));

x_EU(1)= x0;

x_RK = zeros(1,length(t_range));

x_RK(1)= x0;

for k = 1:length(t_range) - 1

x_EU(k+1) = euler_scheme(x_EU(k),dt);

x_RK(k+1) = RK4_scheme(x_RK(k),dt);

end

figure;

plot(t_range,x_EU,'g-o','Markersize',5);

hold on;

plot(t_range,x_RK,'r-o','Markersize',5);

hold on;

t_analytical = t0:0.001:tf;

x_analytical = (exp(t_analytical+log(1/3)))./(1-(exp(t_analytical+log(1/3))));

plot(t_analytical,x_analytical,'k-','linewidth',1);

legend('Euler method','RK4 method','analytical method')

function dxdt = fx(x)

dxdt = (x^2)+x;

end

function x1 = euler_scheme(x0,h)

x1 = x0 + h*fx(x0);

end

function x1 = RK4_scheme(x0,h)

k1 = fx(x0);

k2 = fx(x0 + 0.5*h*k1);

k3 = fx(x0 + 0.5*h*k2);

k4 = fx(x0 + h*k3);

x1 = x0+h*((k1/6)+((k2+k3)/3)+ (k4/6));

end

四阶龙格库塔方法：

%% Question 1 part(f)

clear all;

clc;

tmin = 0;

h = 0.1;

tmax = 10;

xs = [1 2 3];

ys = [0 0 0];

for i = 1:length(xs)

[~,xr(:,2*i-1),yr(:,2*i-1)] = euler_scheme(xs(i),ys(i),h,tmin,tmax);

[~,xr(:,2*i),yr(:,2*i)] = RK4_scheme(xs(i),ys(i),h,tmin,tmax);

end

% [xr,yr] = my_streamline_function(xs, ys, tmin, tmax, h);

figure

legend_context = [];

for i = 1:length(xs) % euler

% plot(xr(1,2*i-1),yr(1,2*i-1))

hold on

plot(xr(:,2*i-1),yr(:,2*i-1),'linewidth',1.5)

xlabel('x')

ylabel('y')

legend_context = [legend_context;['x0=',num2str(xs(i)),', y0=',num2str(ys(i))]];

end

legend(legend_context)

title('Euler Method')

hold off

figure

legend_context2 = [];

for i = 1:length(xs) % RK4

% plot(xr(1,2*i),yr(1,2*i))

hold on

plot(xr(:,2*i),yr(:,2*i),'linewidth',1.5)

xlabel('x')

ylabel('y')

legend_context2 = [legend_context2;['x0=',num2str(xs(i)),', y0=',num2str(ys(i))]];

end

legend(legend_context2)

title('RK4 Method')

hold off

%%

x0 = 0;

y0 = 1;

t = tmin:h:tmax;

for i = 1:length(xs)

[~,x(:,2*i-1),y(:,2*i-1)] = euler_scheme(x0,y0,h,tmin,tmax);

[~,x(:,2*i),y(:,2*i)] = RK4_scheme(x0,y0,h,tmin,tmax);

end

% [x,y] = my_streamline_function(x0, y0, tmin, tmax, h);

t_input = [1,2,3,4,5,6]; % row vector

x_output = zeros(length(t_input),3); % first col:euler, second col:RK4, third col:analytical

y_output = zeros(length(t_input),3);

x_output(:,3) = sin(t_input)';

y_output(:,3) = cos(t_input)';

for k = 1:length(t_input)

if t_input(k) < tmin || t_input(k) > tmax

disp('Error!')

x_output(k,1:2) = NaN;

y_output(k,1:2) = NaN;

else

pos = find(t >= t_input(k));

ind1 = pos(1);

ind2 = pos(1) - 1;

% euler

x_output(k,1) = x(ind2,1)+(x(ind1,1)-x(ind2,1))/(t(ind1)-t(ind2))*(t_input(k)-t(ind2));

y_output(k,1) = y(ind2,1)+(y(ind1,1)-y(ind2,1))/(t(ind1)-t(ind2))*(t_input(k)-t(ind2));

% RK4

x_output(k,2) = x(ind2,2)+(x(ind1,2)-x(ind2,2))/(t(ind1)-t(ind2))*(t_input(k)-t(ind2));

y_output(k,2) = y(ind2,2)+(y(ind1,2)-y(ind2,1))/(t(ind1)-t(ind2))*(t_input(k)-t(ind2));

end

end

x_output

y_output

%% Total function

function [xr,yr] = my_streamline_function(xs, ys, tmin, tmax, h)

end

%% Euler

function [t,x,y] = euler_scheme(x0,y0,dt,t0,tf)

t = t0:dt:tf;

x = zeros(length(t),1);

y = zeros(length(t),1);

x(1) = x0;

y(1) = y0;

for k = 1:length(t) - 1

x(k+1) = x(k) + dt*f(x(k),y(k));

y(k+1) = y(k) + dt*g(x(k),y(k));

end

end

%% RK4

function [t,x,y] = RK4_scheme(x0,y0,h,tmin,tmax)

t = tmin:h:tmax;

x = zeros(length(t),1);

y = zeros(length(t),1);

x(1) = x0;

y(1) = y0;

for k = 1:length(t) - 1

k11 = f(x(k),y(k));

k12 = g(x(k),y(k));

k21 = f(x(k)+h*k11/2,y(k)+h*k12/2);

k22 = g(x(k)+h*k11/2,y(k)+h*k12/2);

k31 = f(x(k)+h*k21/2,y(k)+h*k22/2);

k32 = g(x(k)+h*k21/2,y(k)+h*k22/2);

k41 = f(x(k)+h*k31,y(k)+h*k32);

k42 = g(x(k)+h*k31,y(k)+h*k32);

x(k+1) = x(k) + (k11/6+(k21+k31)/3+k41/6)*h;

y(k+1) = y(k) + (k12/6+(k22+k32)/3+k42/6)*h;

end

end

%% function

function value = f(x,y)

value = y;

end

function value = g(x,y)

value = -x;

end

MATLAB实战——微分方程组的解法之欧拉法与4阶龙格库塔法相关教程

完成使用c++调用matlab生成的.dll库函数的小demo

完成使用c++调用matlab生成的.dll库函数的小demo 完成使用c++调用matlab生成的.dll库函数的小demo 参考了很多人的blog C/C++下调用matlab函数操作说明 C/C++ VS中调用matlab函数的方法 自己动手 按照别人的blog，配置好matlab里边c++编译器配置，我自己用的是

同步异步实战

同步异步实战 源代码下载点击 1.MybatisPlus MyBatis-Plus(简称 MP)是一个 MyBatis 的增强工具，在 MyBatis 的基础上只做增强不做改变，为简化开发、提高效率而生。 2.ORM思想 对象关系映射 (英语：Object Relational Mapping，简称ORM，或O/RM，或O/R map

matlab | 生成低通滤波器

matlab | 生成低通滤波器 引言 生成的低通滤波器如下：(在菲涅尔衍射积分中，我们可认为它是一个衍射面上的衍射孔) 算法 以下算法可直接生成如上图所示的低通滤波器。 r=512,c=r; a=zeros(r,c); a(r/2-r/4:r/2+r/4,c/2-c/4:c/2+c/4)=1; 分析 上述算法难点在于

MATLAB：修改界面左上角Logo图标

MATLAB：修改界面左上角Logo图标 MATLAB 代码 function setLogo(fig, iconPath)% 更改 UI 界面的 Logo ：% setLogo(fig, iconPath)% fig 为窗口句柄% iconPath 为图标存放的路径，类型为 '*.png', '*.jpg', ‘.ico’等 % 参数合法性判断 if ~isvalid(fig), re

如何用函数框架快速开发大型 Web 应用 | 实战

如何用函数框架快速开发大型 Web 应用 | 实战 这是上个月前端早早聊第六届 Serverless 专场的分享。整个分享演示了三个示例，介绍了 Midway Serverless 体系的不同功能，欢迎尝试。 PPT 下载和分享：/midwayjs/midway/tree/resource 大家好

一种快速灰度校正算法(处理亮度不均等情况)(含MATLAB代码)

一种快速灰度校正算法(处理亮度不均等情况)(含MATLAB代码) 文章目录 前言 一、MATLAB代码 二、结果示例 总结 前言 方法来源：[1]高建贞,任明武,杨静宇.一种快速实用的灰度校正算法[J].中国图象图形学报，2002(06):30-34. MATLAB代码：经MATLAB Ra实现

vue2.0 + element-ui 实战项目-axios请求数据(三)

vue2.0 + element-ui 实战项目-axios请求数据(三) 1：进入项目，npm安装 npm install axios --save 2.在main.js下引用axios import axios from 'axios' 3:准备json数据 自己写了一个json数据，放在服务器上，现在要通过vue项目调用数据 http://47.xxx.xx.78

Django实战: Python爬虫爬取链家上海二手房信息，存入数据库并在

Django实战: Python爬虫爬取链家上海二手房信息，存入数据库并在前端显示 今天就带你把它与Python爬虫结合做出个有趣的东西吧。我们将开发这样一个应用，前端用户可以根据行政区划，房厅数和价格区间选择需要爬取的二手房房源信息，后台Python开始爬取数据。