Why is my data is incorrect after looping?

1 visualizzazione (ultimi 30 giorni)
Hello I'm working on this code that solves a set of differential equations and plots the solution and I'm varying the parameter 'k' from 1 to 4 to see how it changes the solution of the problem. For some reason after solving the differential equation inside the loop the u and v data is wrong. To further show this here is what the 'u' data looks like when run for k = 3:
And here is what the data looks like for k = 1:1:4:
As you can see the k=3 column has different data but it should be exactly the same. Can anyone figure out what is going on?
Here is the script:
%Clear command window and workspace
clear
close all
clc
% Fitzhugh-Nagoma model parameters
e=0.01;
k=1:1:4; a=0.1;
i = 0.001;
figure(1);
hold on
u=zeros(500000,4);
v=zeros(500000,4);
t=zeros(100000,1);
% Initial conditions:
u(1)=0.6;
v(1)=0.0;
t(1)=0;
dt=0.001;
%==========================================================================
% Forvard Euler Method, for soluing the ODE
%==========================================================================
for i=1:1:500000
t(i+1)=t(i)+dt;
% u(i+1) = u(i)+ dt*((1/e)*((k(:)*u(i)*(u(i)-a)*(1-u(i)))-v(i)));
u(i+1,:) = u(i,:)+ dt.*((1/e).*((k.*(u(i,:).*(u(i,:)-a).*(1-u(i,:))))-v(i)));
v(i+1,:) = v(i,:)+ dt.*(u(i)-v(i));
end
% Getting the phase plot
figure(1);
upts=(-2:.05:2);
% unullpts=(k(1).*upts.*(upts-a).*(1-upts));
% vnullpts=upts;
for i = 1:numel(k)
unullpts=(k(i).*upts.*(upts-a).*(1-upts));
vnullpts=upts;
plot(upts,unullpts,'black',upts,vnullpts,'black');
hold on
end

Risposta accettata

Torsten
Torsten il 6 Nov 2018
You must set the initial condition for u and v not only for u(1,1) and v(1,1), but for u(1,1:4) and v(1,1:4). And change all u(i) to u(i,:) and v(i) to v(i,:).
I think it is advisable to use a double loop (one from 1 to 4 and the other from 1 to 50000).
Best wishes
Torsten.

Più risposte (0)

Categorie

Scopri di più su Loops and Conditional Statements in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by