[Show all top banners]
Back to: Kurakani General Refresh page to view new replies
 Matlab Help
[VIEWED 3804 TIMES]
SAVE! for ease of future access.
Posted on 09-03-11 5:21 AM     Reply [Subscribe]
Login in to Rate this Post:     0       ?    
 

Hi Matlab Pros,

I need help with this code. I am using RK4 method to solve 2 non-linear equations. yprime is defined quadratically so a, b and c has to be defined. solve these ODEs store the values of y(x) and Ar(x). Equation of Theta is given, plots of x0 vs theta and Ar vs theta should be presented.

function [Theta,xo]=GM_CC
r=0.2;
alpha=50;
x_f=0.209; %inital value
x_0=0.001; %targeted value
n = 209;



Theta =zeros(n,1);
x0=zeros(n,1);

for count=1:1:n
x0(count)=count*((x0_max-x0_min)/n);
Theta(count)=GM_CC(r,x0(count),alpha);

count=count+1;
end
plot(Theta,x0,'m')



function Theta = GM_CC(r,x0,alpha)

% This function simultaneously solves equations A1.10(gives y(x)) and A1.16(gives Ar(x) using 4th order Runge
% Kutta for Counter-Current Outputs of this function include theta (cut-fraction) y, and Ar where Ar is defined as
% Am/Lo
% x0 = mol% O2 in non-permeate stream,x0=(0.209 to 0.001)



h=10^-4; % discretizes space into fine subdivisions for better approximation
xf=0.209; %standard value for oxygen content in air
x=x0+0.1*h; % At time=0, both numerator and denominator are 0, so this offset is added
r=0.2; % ratio of low side to high side pressure.
Am=5.4732; %Calculated Area of Membrane from Halvorson handout (in m^2)

% Defined variables note that the equations include yprime which is solved
% quadratically so we define a,b,c and then the yprime.
a=1-alpha;
b=-1+alpha+(1/r)+(x/r)*(alpha-1);
c=-alpha*(x/r);
yp=(-b+sqrt(b^2-4*a*c))/(2*a);
Ar=0; % Initial value
y = yp; % Stores y' and Ar as one variable
nsteps=ceil((xf-x0)/h); % Defines steps
X(1,1)=x0; Y(1,1)=y'; AR(1,1)=Ar; % Storage vectors for X,Y, and AR

% For loop to solve equations
for i=1:nsteps;
% This ensures we reach xf on last step.
if i==nsteps
h=xf-x;
end

f=rhs1(x,y,x0,r,alpha);
K1=h*f;
f1=rhs1(x+h/2,y+K1/2,x0,r,alpha);
K2=h*f1;
f2=rhs1(x+h/2,y+K2/2,x0,r,alpha);
K3=h*f2;
f3=rhs1(x+h,y+K3,x0,r,alpha);
K4=h*f3;

F=rhs2(x,y,x0,Ar,r,alpha);
J1=h*F;
F1=rhs2(x+h/2,y,x0,Ar+K1/2,r,alpha);
J2=h*F1;
F2=rhs2(x+h/2,y,x0,Ar+K2/2,r,alpha);
J3=h*F2;
F3=rhs2(x+h,y,x0,Ar+K3,r,alpha);
J4=h*F3;

y=y+(K1+2*K2+2*K3+K4)/6;
Ar=Ar+(J1+2*J2+2*J3+J4)/6;

x=x+h;

X(i+1,1)=x;
Y(i+1,1)=y;
AR(i+1,1)=Ar;
end
plot(x0,Ar(i,1),'k')
Theta = (xf-x0)/(Y(i+1,:)-x0);
%Calculates the cut fraction (theta) from final value of y
%display(Theta) %displays calculated cut fraction
%display(Lo) % displays calculated Lo value
%display(Ar) % displays calculated Ar value
%display(y) %displays calculated y value



% The following functions store values that are solved for in the for loop
% shown above
function f=rhs1(x,y,x0,r,alpha) %This function is used to solve for y using equation A1.16
if x==x0
x=x0+0.1*10^-4;
end
% Variables a, b, c, and yp change with time and must be included in this
% function.
r=0.2;
a=1-alpha;
b=-1+alpha+(1/r)+((x/r)*(alpha-1));
c=(-alpha*x)/r;
yp=(-b+sqrt(b^2-4*a*c))/(2*a);

f=((y-x0)*((1-y)*alpha*(x-r*yp)-y*((1-x)-r*(1-yp))))/((x-x0)*((1-x)*alpha*(x-r*yp)-x*((1-x)-r*(1-yp)))); %This is equation A1.16


function F=rhs2(x,y,x0,~,r,alpha) %This funciton is used to solve for Ar using equation A1.10
% Variables a, b, c, and yp change with time and must be included in this
% function.
r=0.2;
a=1-alpha;
b=-1+alpha+(1/r)+((x/r)*(alpha-1));
c=(-alpha*x)/r;
yp2=(-b+sqrt(b^2-4*a*c))/(2*a);

F=(-((y-x0)/(x-y)))/((((1-x)*alpha*(x-r*yp2)-x*((1-x)-r*(1-yp2))))); %This is equation A1.10

Please help me or give suggestions for the code.

 


Please Log in! to be able to reply! If you don't have a login, please register here.

YOU CAN ALSO



IN ORDER TO POST!




Within last 365 days
Recommended Popular Threads Controvertial Threads
शीर्षक जे पनि हुन सक्छ।
TPS Re-registration case still pending ..
tesla stock OMG !!
What are your first memories of when Nepal Television Began?
TPS for Nepal likely to extend next week
Anybody gotten the TPS EAD extension alert notice (i797) thing? online or via post?
MAGA #FAFO
Embassy of Nepal might be able to help extend TPS for Nepal
I hope all the fake Nepali refugee get deported
Basnet or Basnyat ??
TPS EAD auto extended to June 2025 or just TPS?
Those who are in TPS, what’s your backup plan?
nrn citizenship
Toilet paper or water?
Sajha has turned into MAGATs nest
ChatSansar.com Naya Nepal Chat
Nas and The Bokas: Coming to a Night Club near you
ढ्याउ गर्दा दसैँको खसी गनाउच
Mamta kafle bhatt is still missing
Who is hottest nepali female?
Nas and The Bokas: Coming to a Night Club near you
Looking for girl
Who is hottest nepali female?
NOTE: The opinions here represent the opinions of the individual posters, and not of Sajha.com. It is not possible for sajha.com to monitor all the postings, since sajha.com merely seeks to provide a cyber location for discussing ideas and concerns related to Nepal and the Nepalis. Please send an email to admin@sajha.com using a valid email address if you want any posting to be considered for deletion. Your request will be handled on a one to one basis. Sajha.com is a service please don't abuse it. - Thanks.

Sajha.com Privacy Policy

Like us in Facebook!

↑ Back to Top
free counters