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

3 visualizaciones (últimos 30 días)
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

Respuesta aceptada

William Rose
William Rose el 3 de Mayo de 2021
Editada: William Rose el 3 de Mayo de 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 comentarios
William Rose
William Rose el 3 de Mayo de 2021
@Anil Maddipatla, Thank you for posting some code .
Please answer my original question: What is the equation for the derivative of the funciton you want to compute?
You wrote in the original post:
Is c(t) the derivative of y(t), or is c(t) the solution? I suspect it is the solution. If it is the solution, then what is the derivative, ? I gave a suggestion in my previous answer.
You said you do not fully understand the code you posted, and you said you got the code from somewhere else, and edited it. I too do not fully understand it. The top part of your posting is a function, Fn_Euler(), which returns vectors x and y. Fn_Euler() will not work. For one thing, the function should be passed as a handle, with an @ sign. Also, the function f is defined inside the for loop , after it is used the first time, so the function f that is passed (if it were passed corectly, which it is not) will only be used once, on the first time through the foor loop. Also, the parameters w, v, and lambda are not passed to the function, and are not defined inside the function.
The second part of your posting appears to be a main program. It does not call Fn_Euler. This script attemts to do the integration by the Euler method. It will not work as desired, partly because f is not defined as a function. It is simply a variable to which a value is assigned. After the value is assigned, it is not updated.
The script defines the value of W at every time point (1 year intervals). The values change abruptly at certain times. Are you sure the vales of w really change in this manner?
Please answer my first question first: is c(t) the solution, or is c(t) the derivative of y(t)?
Anil Maddipatla
Anil Maddipatla el 3 de Mayo de 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

Iniciar sesión para comentar.

Más respuestas (1)

William Rose
William Rose el 3 de Mayo de 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.

Categorías

Más información sobre Programming en Help Center y File Exchange.

Productos


Versión

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by