MATLAB Answers

How to code for this equation based on euler's, Huen's and 4th order Range Kutta method of ODE on MATLAB. Please Help.

12 views (last 30 days)
Hi,
I have to code for computing & plotting this equation in Euler's Method, Huen's Method & 4th Order Range Kutta Method separately.
There are various values of loading "w" in gram per year and corresponding values of "t" in years. "Co" value is given. Lambda shall be calculated by the equation lambda = ((Q/V)+k+(u/H)), where the values of Q, V, u, H, k are available. Time step is given as 1. I am very new to MATLAB and I need to submit the plot and output. Please help.
Regards,
Anil

Accepted Answer

William Rose
William Rose on 3 May 2021
Edited: William Rose on 3 May 2021
What have you tried so far, and what errors do you get?
The three methods you have cited are for solving a differential equation. However, your question does not show a differential equation.
It's like Differential Equation Jeopardy: you have provided the answer, but what is the question?
It looks to me like the question is: Solve the differential equation
,
with the initial condition .
Is that correct? Please specify the time period for which you must solve it (start at t=0, end at what time?), and the values of and all the other constants you mentioned.
Once you have the differential equation, write a script to implement forward Euler method, and the Heun method, and then the ode45() approach.
For all three methods, you must create a function that returns the value of the derivative, given the value of c and the time.
For the ode45() method, you will pass the name of the derivative-returning function to ode45(), and ode45() will then figure out the solution. ode45() also needs to know the start time, the end time, and the initial condition. You do not need to specify a step size, because ode45() will figure out a good step size as it goes, and will adapt the step size as it goes, to maintain a good balance between accuracy and efficiency.
For the other two methods, you will have to specify a step size, and you will write a for loop that goes step by step across the full time interval. You will call the function that computes the derivative once per step, for Euler. You will call it twice per step, for the Heun method.
That should be enough to get you going. Write some code, and if it does not work, post what you have written, and any error messages you get. Good luck.
  3 Comments
Anil Maddipatla
Anil Maddipatla on 3 May 2021
Hi Rose,
Thank you. Yes. You are right. C(t) is the solution. All of which you mentioned in your reply is exactly right. After your reply, I have tried the code in a simpler way, which luckily worked. I have posted the code below:
For Euler Method: t =[0:60]; %Time in Years%
load('W_Data.mat') %loading in g per year% V = 2900000000; %Volume in m3% v= 12; %settling velocity in m per year% Q = 1250000000; %values of Outflow in m3 per year% H = 33; %values of mean depth in m% h=1; %values of Time Step in years% k=0; %values of Rate of Decay per year% Lambda = (Q/V)+k+(v/H); %Eigen Value lambda per year% c(1)= 0.0174; %the value of intial concentration in milligram per year% for i=1:1:60 c(i+1)= c(i)+((w(1,i)/V)-(Lambda*c(i)))*1 end t=[1930:1990] plot(t,c,'-o') title(' Euler Method')
For Huen Method: t = [0:60]; %Time in year% load('W_Data.mat') %loading in g per year% V = 2900000000; %Volume in m3% v = 12; %settling velocity in m per year% Q = 1250000000; %values of Outflow in m3 per year% H = 33; %values of mean depth in m% h = 1; %values of time step size in years% k=0; %values of Rate of Decay per year% Lambda = (Q/V)+k+(v/H); %Eigen Value lambda per year% c(1) = 0.0174; %the value of intial concentration in milligram per year% for i=1:1:60 c0(i+1) = c(i)+((w(1,i)/V)-(Lambda*c(i)))*h; c(i+1) = c(i)+(((w(1,i)/V)-(Lambda*c(i)))+((w(1,i)/V)-(Lambda*c0(i+1))))/2 end t=[1930:1990]; plot(t,c,'-o') title(' Huen Method')
For Runge Kutta 4th derivative method: t = [0:60]; %time in year% load('W_Data.mat') %loading in g per year% V = 2900000000; %Volume in m3% v = 12; %settling velocity in m per year% Q = 1250000000; %values of Outflow in m3 per year% H = 33; %values of mean depth in m% h = 1; %values of time step size in years% k=0; %values of Rate of Decay per year% Lambda = (Q/V)+k+(v/H); %Eigen Value lambda per year% c(1) = 0.0174; %the value of intial concentration in milligram per year% for i=1:1:60 k1(i)= (w(1,i)/V)-(Lambda*c(i)); k2(i)= (w(1,i)/V)-(Lambda*c(i)+(k1(i)/2)); k3(i)= (w(1,i)/V)-(Lambda*c(i)+(k2(i)/2)); k4(i)= (w(1,i)/V)-(Lambda*c(i)+(k3(i))); c(i+1)= c(i)+ (((k1(i))+(2*k2(i))+(2*k3(i))+(k4(i)))/6) end t=[1930:1990]; plot(t,c,'-o') title('RK4 Method')
I have generated individual plots for these methods. But, if I want to generate all the graphs on a combined plot, can you guide me how to do it please.
Once again thanks for your replies. They helped a lot.
Regards,
Anil

Sign in to comment.

More Answers (1)

William Rose
William Rose on 3 May 2021
You're welcome.
To plot all three solutions on one plot, I recommend that you name the solutions c1, c2, c3 (for Euler, Heun, and R-K). You can plot them on the same plot as follows:
figure;
plot(t,c1,'rx-',t,c2,'go-',t,c3,'b+-');
legend('Euler','Heun','R-K');
xlabel('Time (years)');
ylabel('Concentration');
title('Concentration versus Time: Three Solutions');
I think the formula for the derivative which I posted earlier was not quite correct. My new estimate for the derivative is
The derivative above is consistent with your Euler update formula:
c(i+1)= c(i)+((w(1,i)/V)-(Lambda*c(i)))*1
which implies that
dc/dt=w(1,i)/V-Lambda*c(i),
and you have used h=1.

Community Treasure Hunt

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

Start Hunting!