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开始爬取数据。