Some one check this code of solving ODEs . Is it correct or not?

5 vues (au cours des 30 derniers jours)
Akhtar Jan
Akhtar Jan le 20 Nov 2022
Modifié(e) : Akhtar Jan le 21 Nov 2022
n = 500;
tmin = 0;
tmax = 8;
h = tmax/n;
t = linspace(tmin, tmax, n)
E = 1;
x = zeros(1, n);
x(1) = 8;
y = zeros(1, n);
y(1) = 0;
z = zeros(1, n);
z(1) = 0;
w = zeros(1, n);
w(1) = 0;
f1 =@(t, x,y , z, w) -k1*x*E+k_1*y-k3*x*y+k_3*z;
f2 =@(t, x, y, z, w) k1*x*E-(k_1+k2)*y-k3*x*y*+(k_3+k4)*z;
K1 = f1(t(i), x(i), y(i), z(i), w(i));
l1 = f2(t(i), x(i), y(i), z(i), w(i));
K2 = f1(t(i)+h/2, x(i)+K1*h/2, y(i)+l1*h/2, z(i)+m1*h/2, w(i)+n1*h/2);
l2 = f2(t(i)+h/2, x(i)+K1*h/2, y(i)+l1*h/2, z(i)+m1*h/2, w(i)+n1*h/2);
K3 = f1(t(i)+h/2, x(i)+K2*h/2, y(i)+l2*h/2, z(i)+m2*h/2, w(i)+n2*h/2);
l3 = f2(t(i)+h/2, x(i)+K2*h/2, y(i)+l2*h/2, z(i)+m2*h/2, w(i)+n2*h/2);
K4 = f1(t(i)+h, x(i)+K3*h, y(i)+l3*h, z(i)+m3*h, w(i)+n3*h);
l4 = f2(t(i)+h, x(i)+K3*h, y(i)+l3*h, z(i)+m3*h, w(i)+n3*h);
x(i+1) = x(i)+h/6*(K1 +2*K2 + 2*K3 +K4);
y(i+1) = y(i)+h/6*(l1 +2*l2 + 2*l3 +l4);
end
for i= 1:4
t0 = tmin+i*h;
part11 =55.0*f1(t(4), x(4), y(4), z(4), w(4))-59.0*f1(t(3), x(3), y(3), z(3), w(3))+37.0*f1(t(2), x(2), y(2), z(2), w(2))-9.0*f1(t(1), x(1), y(1), z(1), w(1));
x011 = x(4)+h*(part11)/24.0;
part22 =55.0*f2(t(4), x(4), y(4), z(4), w(4))-59.0*f2(t(3), x(3), y(3), z(3), w(3))+37.0*f2(t(2), x(2), y(2), z(2), w(2))-9.0*f2(t(1), x(1), y(1), z(1), w(1));
y011 = y(4)+h*(part22)/24.0;
part33 =55.0*f3(t(4), x(4), y(4), z(4), w(4))-59.0*f3(t(3), x(3), y(3), z(3), w(3))+37.0*f3(t(2), x(2), y(2), z(2), w(2))-9.0*f3(t(1), x(1), y(1), z(1), w(1));
z011 = z(4)+h*(part33)/24.0;
part44 =55.0*f4(t(4), x(4), y(4), z(4), w(4))-59.0*f4(t(3), x(3), y(3), z(3), w(3))+37.0*f4(t(2), x(2), y(2), z(2), w(2))-9.0*f4(t(1), x(1), y(1), z(1), w(1));
w011 = w(4)+h*(part44)/24.0;
part11 = 9.0*f1(t0, x011, y011, z011,w011)+19*f1(t(4), x(4), y(4), z(4))-5.0*f1(t(3), x(3), y(3), z(3))+ f1(t(2), x(2), y(2), z(2));
x021 = x(4) + h*(part11)/24;
part22 = 9.0*f2(t0, x011, y011, z011,w011)+19*f2(t(4), x(4), y(4), z(4),w(4))-5.0*f2(t(3), x(3), y(3), z(3),w(3))+ f2(t(2), x(2), y(2), z(2), w(2));
y011 = y(4) + h*(part22)/24;
part33 = 9.0*f3(t0, x011, y011, z011,w011)+19*f3(t(4), x(4), y(4), z(4),w(4))-5.0*f3(t(3), x(3), y(3), z(3),w(3))+ f3(t(2), x(2), y(2), z(2),w(2));
z011 = z(4) + h*(part33)/24;
part44 = 9.0*f4(t0, x011, y011, z011,w011)+19*f4(t(4), x(4), y(4), z(4),w(4))-5.0*f4(t(3), x(3), y(3), z(3),w(3))+ f4(t(2), x(2), y(2), z(2), w(2));
w011 = w(4) + h*(part44)/24;
for j = 1:3
t(j) = t(j+1);
x(j)= x(j+1);
y(j)= y(j+1);
z(j)= z(j+1);
w(j)= w(j+1);
end
t(4) = t0;
x(4) = x011;
y(4) = y011;
end
pt1 = plot(t, x,'r');
hold on
pt2 = plot(t, y,'g');
pt3 = plot(t, z , 'b');
pt4 = plot(t, w,'k');
  3 commentaires
Akhtar Jan
Akhtar Jan le 20 Nov 2022
Sir please check out the graph at the start the graph is not a smooth curve
John D'Errico
John D'Errico le 20 Nov 2022
I think you don't understand. Your code is quite difficult to read. Quite lengthy. Not a single comment in it. Lots of stuff hard coded. And that makes it almost necessary to decode what you did on this lengthy mess of code. What for example, do variables like part 33 and part44 do? (I'm just picking two of them as an example.)
On top of that, it is not even obvious what differential equation you are trying to solve, since you have hard coded functions directly in the code.
When you test a code, you need to test that code on something where you ABSOLUTELY KNOW the solution.
  1. Apply the code to a trivial problem. Verify you get the correct solution.
  2. Apply the code to a slightly more difficult problem, but again, one where you know the solution. Then see if reducing the step size properly reduces the error.
Those steps would require that it is easily possible to change the differential equation. Of course, any high quality code would have that be trivial, but this is cearly an assignment by a student, as otherwise one would never write their own ODE solver, when far better tools are already available.
But decoding what you did could easily take someone hours to do. And that makes it far less likely you will get a helpful answer, certainly not before this assignment of yours is due. I would STRONGLY suggest that in the future, you need to learn to write better code. (Sorry, but what can I say?) NEVER hard code your function directly into the solver. Learn to use functions. Learn to pass parameters into those functions, so you can easily change the problem to test it out. Had you done things like that, helping you would be far easier now. OF course then, you probably would have been able to perform the tests I suggested above, and know if your code is working properly anyway.

Connectez-vous pour commenter.

Réponse acceptée

Torsten
Torsten le 20 Nov 2022
Now it's smooth. But you still don't know whether the results are correct if you don't test your code with differential equations you know the solution of.
n = 5000;
tmin = 0;
tmax = 8;
h = tmax/n;
t = linspace(tmin, tmax, n);
k1 = 3;
k_1 = 2;
k2 = 7;
k3 = 1;
k_3 = .1;
k4 = 2.5;
E = 1;
x = zeros(1, n);
x(1) = 8;
y = zeros(1, n);
y(1) = 0;
z = zeros(1, n);
z(1) = 0;
w = zeros(1, n);
w(1) = 0;
f1 =@(t, x,y , z, w) -k1*x*E+k_1*y-k3*x*y+k_3*z;
f2 =@(t, x, y, z, w) k1*x*E-(k_1+k2)*y-k3*x*y*+(k_3+k4)*z;
f3 =@(t, x, y, z, w) k3*x*y-(k_3+k4)*z;
f4 =@(t, x, y, z, w) k2*y+k4*z;
for i=1:n-1
K1 = f1(t(i), x(i), y(i), z(i), w(i));
l1 = f2(t(i), x(i), y(i), z(i), w(i));
m1 = f3(t(i), x(i), y(i), z(i), w(i));
n1 = f4(t(i), x(i), y(i), z(i), w(i));
K2 = f1(t(i)+h/2, x(i)+K1*h/2, y(i)+l1*h/2, z(i)+m1*h/2, w(i)+n1*h/2);
l2 = f2(t(i)+h/2, x(i)+K1*h/2, y(i)+l1*h/2, z(i)+m1*h/2, w(i)+n1*h/2);
m2 = f3(t(i)+h/2, x(i)+K1*h/2, y(i)+l1*h/2, z(i)+m1*h/2, w(i)+n1*h/2);
n2 = f4(t(i)+h/2, x(i)+K1*h/2, y(i)+l1*h/2, z(i)+m1*h/2, w(i)+n1*h/2);
K3 = f1(t(i)+h/2, x(i)+K2*h/2, y(i)+l2*h/2, z(i)+m2*h/2, w(i)+n2*h/2);
l3 = f2(t(i)+h/2, x(i)+K2*h/2, y(i)+l2*h/2, z(i)+m2*h/2, w(i)+n2*h/2);
m3 = f3(t(i)+h/2, x(i)+K2*h/2, y(i)+l2*h/2, z(i)+m2*h/2, w(i)+n2*h/2);
n3 = f4(t(i)+h/2, x(i)+K2*h/2, y(i)+l2*h/2, z(i)+m2*h/2, w(i)+n2*h/2);
K4 = f1(t(i)+h, x(i)+K3*h, y(i)+l3*h, z(i)+m3*h, w(i)+n3*h);
l4 = f2(t(i)+h, x(i)+K3*h, y(i)+l3*h, z(i)+m3*h, w(i)+n3*h);
m4 = f3(t(i)+h, x(i)+K3*h, y(i)+l3*h, z(i)+m3*h, w(i)+n3*h);
n4 = f4(t(i)+h, x(i)+K3*h, y(i)+l3*h, z(i)+m3*h, w(i)+n3*h);
x(i+1) = x(i)+h/6*(K1 +2*K2 + 2*K3 +K4);
y(i+1) = y(i)+h/6*(l1 +2*l2 + 2*l3 +l4);
z(i+1) = z(i)+h/6*(m1 +2*m2 + 2*m3 +m4);
w(i+1) = w(i)+h/6*(n1 +2*n2 + 2*n3 +n4);
end
for i= 1:4
t0 = tmin+i*h;
part11 =55.0*f1(t(4), x(4), y(4), z(4), w(4))-59.0*f1(t(3), x(3), y(3), z(3), w(3))+37.0*f1(t(2), x(2), y(2), z(2), w(2))-9.0*f1(t(1), x(1), y(1), z(1), w(1));
x011 = x(4)+h*(part11)/24.0;
part22 =55.0*f2(t(4), x(4), y(4), z(4), w(4))-59.0*f2(t(3), x(3), y(3), z(3), w(3))+37.0*f2(t(2), x(2), y(2), z(2), w(2))-9.0*f2(t(1), x(1), y(1), z(1), w(1));
y011 = y(4)+h*(part22)/24.0;
part33 =55.0*f3(t(4), x(4), y(4), z(4), w(4))-59.0*f3(t(3), x(3), y(3), z(3), w(3))+37.0*f3(t(2), x(2), y(2), z(2), w(2))-9.0*f3(t(1), x(1), y(1), z(1), w(1));
z011 = z(4)+h*(part33)/24.0;
part44 =55.0*f4(t(4), x(4), y(4), z(4), w(4))-59.0*f4(t(3), x(3), y(3), z(3), w(3))+37.0*f4(t(2), x(2), y(2), z(2), w(2))-9.0*f4(t(1), x(1), y(1), z(1), w(1));
w011 = w(4)+h*(part44)/24.0;
part11 = 9.0*f1(t0, x011, y011, z011,w011)+19*f1(t(4), x(4), y(4), z(4))-5.0*f1(t(3), x(3), y(3), z(3))+ f1(t(2), x(2), y(2), z(2));
x021 = x(4) + h*(part11)/24;
part22 = 9.0*f2(t0, x011, y011, z011,w011)+19*f2(t(4), x(4), y(4), z(4),w(4))-5.0*f2(t(3), x(3), y(3), z(3),w(3))+ f2(t(2), x(2), y(2), z(2), w(2));
y011 = y(4) + h*(part22)/24;
part33 = 9.0*f3(t0, x011, y011, z011,w011)+19*f3(t(4), x(4), y(4), z(4),w(4))-5.0*f3(t(3), x(3), y(3), z(3),w(3))+ f3(t(2), x(2), y(2), z(2),w(2));
z011 = z(4) + h*(part33)/24;
part44 = 9.0*f4(t0, x011, y011, z011,w011)+19*f4(t(4), x(4), y(4), z(4),w(4))-5.0*f4(t(3), x(3), y(3), z(3),w(3))+ f4(t(2), x(2), y(2), z(2), w(2));
w011 = w(4) + h*(part44)/24;
for j = 1:3
t(j) = t(j+1);
x(j)= x(j+1);
y(j)= y(j+1);
z(j)= z(j+1);
w(j)= w(j+1);
end
t(4) = t0;
x(4) = x011;
y(4) = y011;
z(4) = z011;
w(4) = w011;
end
pt1 = plot(t, x,'r');
hold on
pt2 = plot(t, y,'g');
pt3 = plot(t, z , 'b');
pt4 = plot(t, w,'k');
  1 commentaire
Akhtar Jan
Akhtar Jan le 20 Nov 2022
Thank you sir you make it possible i have checked my set of ODEs our solution is correct now Thanks🙂

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Tags

Produits


Version

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by