How to replace ode45 with Eulers method?

8 次查看（过去 30 天）
Westin Messer
Westin Messer 2018-10-8
编辑： Torsten 2018-10-9
This code currently works perfectly but I'm being asked to solve the differential equation using Euler's method instead of ode45 and I'm not sure how to do that. Could anyone help me replace ode45 with Euler's method? Thanks!
clear
close all
clc
tic
% INPUT SECTION ========================================================
% External current stimulus [Iext = 1.0]
Iext = 0.33;
% Initial conditions: y0 = [ v(0) = -2.8, w(0) = -1.8]
y0 = [-2.8; -1.8];
% Simulation time [tMin = 0 tMax = 200]
tSpan = [0 200];
% Phase space plot dimensions: X axis [-3 3] / Y axis [-2 3]
Lx = [-3 3]; Ly = [-2 3];
% Phase plot ticks
ticksX = Lx(1):1:Lx(2);
ticksY = Ly(1):1:Ly(2);
% Phase space setup for vector field: number of vectors nX [20]
nX = 20;
% Fitzhugh-Nagoma model parameters
a = 0.7; b = 0.8; c = 12.5;
% CALCULATION SECTION ===================================================
% K use to pass variables to function FNode
K = [Iext; a; b; c];
% Solve differential equations using ode45
[t,y] = ode45(@(t,y) FNode(t,y,K), tSpan,y0);
% Vector field
x1 = linspace(Lx(1),Lx(2),nX);
x2 = linspace(Ly(1),Ly(2),nX);
[xx, yy] = meshgrid(x1,x2);
f = (xx - xx.^3/3 - yy + Iext);
g = (1/c).*(xx + a - b.*yy);
fs = f./sqrt(f.^2 + g.^2); % unit vectors
gs = g./sqrt(f.^2 +g.^2);
% Critical point - output to Command Window
syms p
Sp = vpasolve(p-p^3/3-(p+a)/b + Iext == 0,p,[-3 3]);
Sq = (Sp+a)/b;
Sp = double(Sp); Sq = double(Sq);
disp('Critical point');
fprintf(' v_C = %2.2f\n', Sp);
disp(' ')
fprintf(' w_C = %2.2f\n', Sq);
% GRAPHICS SECTION=======================================================
FS = 14; % fontsize
% Phase space plot: v vs w ------------------------------------------
figure(1)
pos = [0.35 0.05 0.29 0.39];
set(gcf,'Units','normalized');
set(gcf,'Position',pos);
set(gcf,'color','w');
hold on
box on
% VECTOR FIELD
hq = quiver(xx,yy,fs,gs);
set(hq,'color',[0.2 0.2 0.2],'AutoScaleFactor',0.6);
set(gca,'fontsize',FS)
xlim(Lx)
ylim(Ly)
set(gca,'xtick',ticksX);
set(gca,'ytick',ticksY);
grid on
xlabel('membrane potential v'); ylabel('recovery variable w');
% v nullcline
v = linspace(Lx(1),Lx(2),200);
w = (v - v.^3/3 + Iext);
xP = v; yP = w;
plot(xP,yP,'r','linewidth',1.5)
% r nullcline
w = (v + a)/b;
xP = v; yP = w;
plot(xP,yP,'m','linewidth',1.5)
% Phase portrait trajectory
xP = y(:,1); yP = y(:,2);
plot(xP,yP,'b','linewidth',2)
xP = y(1,1); yP = y(1,2); % initial conditions: start of trajectory
Hplot = plot(xP,yP,'o');
set(Hplot,'markersize',8,'markerfacecolor',[0 1 0],'markeredgecolor',[0 1 0])
xP = y(end,1); yP = y(end,2); % end of trajectory
Hplot = plot(xP,yP,'o');
set(Hplot,'markersize',8,'markerfacecolor',[1 0 0],'markeredgecolor',[1 0 0])
tm1 = 'I_{ext} = ';
tm2 = num2str(Iext,'%3.3f');
tm3 = ' v_C = ';
tm4 = num2str(Sp,'%3.2f');
tm5 = ' w_C = ';
tm6 = num2str(Sq,'%3.2f');
tm = [tm1 tm2 tm3 tm4 tm5 tm6];
hT = title(tm,'FontName','Courier');
% Time evolution of v and w
figure(2)
pos = [0.05 0.05 0.29 0.29];
set(gcf,'Units','normalized');
set(gcf,'Position',pos);
set(gcf,'color','w');
xP = t; yP = y(:,1); % v
plot(xP,yP,'b','linewidth',2)
hold on
yP = y(:,2); % w
plot(xP,yP,'r','linewidth',2)
legend('v','w','location','south','orientation','horizontal')
xlabel('t')
ylabel('v & w')
title(tm,'fontName','Courier')
grid on
set(gca,'fontsize',FS)
box on
disp(' ')
toc
% FUNCTIONS ===========================================================
function dydt = FNode(t,y,K)
a = K(2); b = K(3); c = K(4); Iext = K(1);
%Iext = 0.003 * t;
%if t > 50; Iext = 0.5; end
%if t > 55; Iext = 0; end
dydt = [(y(1) - y(1)^3/3 - y(2) + Iext); (1/c)*(y(1) + a - b*y(2))];
end

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2018-10-9
编辑：Torsten 2018-10-9
...
% Solve differential equations using ode45
[t,y] = euler(@FNode,tSpan,y0,K);
...
function [tsol, ysol] = euler(fun, tspan, y0, K)
n = floor(20*tspan(2));
%Compute h from the time grid.
h = (tspan(2) - tspan(1))/n;
%How many equations?
m = length(y0);
%Initilize t and y
t = tspan(1);
y = reshape(y0, 1, m);
% Store initial conditions
tsol = t;
ysol = y;
% Euler loop
for i = 1: n
% Take the Euler step into the temporary variable
y1 = y + h * fun(t, y, K).';
% Update the temporary variables
t = t + h;
y = y1;
% Update the solution vector
tsol = [tsol ; t];
ysol = [ysol ; y];
end
end

更多回答（2 个）

Jan
Jan 2018-10-8
编辑：Jan 2018-10-8
I recommend to ask an internet search engine for "Matlab euler method". You will find e.g.:
Does this help already? Try to implement it. It is not tricky to expand this to a 2D y. Post what you have tried so far and ask a specific question in case of problems.
KSSV
KSSV 2018-10-8
Read about ode23
  16 个评论
显示 14更早的评论
Jan
Jan 2018-10-8
@Westin Messer: You can define the time either by a vector, or by a range and a number of steps. This is fairly equivalent.
Westin Messer
Westin Messer 2018-10-8
This is what I have now but I'm getting a lot of errors. Am I getting close?
clear
close all
clc
tic
% INPUT SECTION ========================================================
% External current stimulus [Iext = 1.0]
Iext = 0.33;
% Initial conditions: y0 = [ v(0) = -2.8, w(0) = -1.8]
y0 = [-2.8; -1.8];
% Simulation time [tMin = 0 tMax = 200]
tSpan = [0 200];
% Phase space plot dimensions: X axis [-3 3] / Y axis [-2 3]
Lx = [-3 3]; Ly = [-2 3];
% Phase plot ticks
ticksX = Lx(1):1:Lx(2);
ticksY = Ly(1):1:Ly(2);
% Phase space setup for vector field: number of vectors nX [20]
nX = 20;
% Fitzhugh-Nagoma model parameters
a = 0.7; b = 0.8; c = 12.5;
% CALCULATION SECTION ===================================================
% K use to pass variables to function FNode
K = [Iext; a; b; c];
% Solve differential equations using ode45
[t,y] = euler(@FNode,tSpan,y0);
% Vector field
x1 = linspace(Lx(1),Lx(2),nX);
x2 = linspace(Ly(1),Ly(2),nX);
[xx, yy] = meshgrid(x1,x2);
f = (xx - xx.^3/3 - yy + Iext);
g = (1/c).*(xx + a - b.*yy);
fs = f./sqrt(f.^2 + g.^2); % unit vectors
gs = g./sqrt(f.^2 +g.^2);
% Critical point - output to Command Window
syms p
Sp = vpasolve(p-p^3/3-(p+a)/b + Iext == 0,p,[-3 3]);
Sq = (Sp+a)/b;
Sp = double(Sp); Sq = double(Sq);
disp('Critical point');
fprintf(' v_C = %2.2f\n', Sp);
disp(' ')
fprintf(' w_C = %2.2f\n', Sq);
% GRAPHICS SECTION=======================================================
FS = 14; % fontsize
% Phase space plot: v vs w ------------------------------------------
figure(1)
pos = [0.35 0.05 0.29 0.39];
set(gcf,'Units','normalized');
set(gcf,'Position',pos);
set(gcf,'color','w');
hold on
box on
% VECTOR FIELD
hq = quiver(xx,yy,fs,gs);
set(hq,'color',[0.2 0.2 0.2],'AutoScaleFactor',0.6);
set(gca,'fontsize',FS)
xlim(Lx)
ylim(Ly)
set(gca,'xtick',ticksX);
set(gca,'ytick',ticksY);
grid on
xlabel('membrane potential v'); ylabel('recovery variable w');
% v nullcline
v = linspace(Lx(1),Lx(2),200);
w = (v - v.^3/3 + Iext);
xP = v; yP = w;
plot(xP,yP,'r','linewidth',1.5)
% r nullcline
w = (v + a)/b;
xP = v; yP = w;
plot(xP,yP,'m','linewidth',1.5)
% Phase portrait trajectory
xP = y(:,1); yP = y(:,2);
plot(xP,yP,'b','linewidth',2)
xP = y(1,1); yP = y(1,2); % initial conditions: start of trajectory
Hplot = plot(xP,yP,'o');
set(Hplot,'markersize',8,'markerfacecolor',[0 1 0],'markeredgecolor',[0 1 0])
xP = y(end,1); yP = y(end,2); % end of trajectory
Hplot = plot(xP,yP,'o');
set(Hplot,'markersize',8,'markerfacecolor',[1 0 0],'markeredgecolor',[1 0 0])
tm1 = 'I_{ext} = ';
tm2 = num2str(Iext,'%3.3f');
tm3 = ' v_C = ';
tm4 = num2str(Sp,'%3.2f');
tm5 = ' w_C = ';
tm6 = num2str(Sq,'%3.2f');
tm = [tm1 tm2 tm3 tm4 tm5 tm6];
hT = title(tm,'FontName','Courier');
% Time evolution of v and w
figure(2)
pos = [0.05 0.05 0.29 0.29];
set(gcf,'Units','normalized');
set(gcf,'Position',pos);
set(gcf,'color','w');
xP = t; yP = y(:,1); % v
plot(xP,yP,'b','linewidth',2)
hold on
yP = y(:,2); % w
plot(xP,yP,'r','linewidth',2)
legend('v','w','location','south','orientation','horizontal')
xlabel('t')
ylabel('v & w')
title(tm,'fontName','Courier')
grid on
set(gca,'fontsize',FS)
box on
disp(' ')
toc
% FUNCTIONS ===========================================================
function dydt = FNode(t,y,K)
a = K(2); b = K(3); c = K(4); Iext = K(1);
%Iext = 0.003 * t;
%if t > 50; Iext = 0.5; end
%if t > 55; Iext = 0; end
dydt = [(y(1) - y(1)^3/3 - y(2) + Iext); (1/c)*(y(1) + a - b*y(2))];
end
function [tsol, ysol] = euler(fun, tspan, y0)
if nargin(fun) ~=2
error('fun must take two inputs, t and y.');
end
if ~all(size(y0) == size(fun(0, y0)))
error('You have not passed appropriate fun or y_0.');
end
%Set up the time grid. ***NOTE THE n+1***
tsol = linspace(0, T, n+1);
%Compute h from the time grid.
h = tsol(2) - tsol(1);
%Orient tgrid as a column vector.
tsol = reshape(tsol, n+1, 1);
%How many equations?
m = length(y0);
%Orient y0 as a row vector.
y0 = reshape(y0, 1, m);
% Preallocate an array to hold the approximate solution. Each row
% corresponds to a point in the time grid.
ysol = zeros(n+1, m);
% Set the initial conditions.
ysol(1,:) = y0;
% Euler loop
for i = 1: n
% Store the point in time as a temporary variable
t_i = tsol(i);
% Take the Euler step into the temporary variable
y_1 = y0 + h * fun(t_i, y0);
% Store the Euler step
ysol(i+1,:) = y_1;
% Update the temporary variable
y0 = y_1;
end
end

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by